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ABSTRACT 

Stellar population synthesis techniques for predicting the observable light emitted 
by a stellar population have extensive applications in numerous areas of astronomy. 
However, accurate predictions for small populations of young stars, such as those found 
in individual star clusters, star-forming dwarf galaxies, and small segments of spiral 
galaxies, require that the population be treated stochastically. Conversely, accurate 
deductions of the properties of such objects also requires consideration of stochasticity. 
Here we describe a comprehensive suite of modular, open-source software tools for 
tackling these related problems. These include: a greatly-enhanced version of the slug 
code introduced by da Silva et al. (2012), which computes spectra and photometry for 
stochastically- or deterministically-sampled stellar populations with nearly-arbitrary 
star formation histories, clustering properties, and initial mass functions; cloudy_slug, 
a tool that automatically couples slug-computed spectra with the cloudy radiative 
transfer code in order to predict stochastic nebular emission; bayesphot, a general- 
purpose tool for performing Bayesian inference on the physical properties of stellar 
systems based on unresolved photometry; and cluster_slug and SFR_slug, a pair 
of tools that use bayesphot on a library of slug models to compute the mass, age, 
and extinction of mono-age star clusters, and the star formation rate of galaxies, 
respectively. The latter two tools make use of an extensive library of pre-computed 
stellar population models, which are included the software. The complete package is 
available at http: //www. slugsps . com. 

Key words: methods: statistical; galaxies: star clusters: general; galaxies: stellar 
content; stars: formation; methods: numerical; techniques: photometric 


1 INTRODUCTION 

Stellar population synthesis (SPS) is a critical tool that al¬ 
lows us to link the observed light we receive from unresolved 
stellar populations to the physical properties (e.g. mass, 
age) of the emitting stars. Reflecting this importance, over 
the years numerous research groups have written and dis¬ 
tributed SPS codes such as starburst99 (Leitherer et al. 
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1999; Vazquez V Leitherer 2005), fsps (Conroy et al. 2009; 
Conroy V Gunn 2010), pegase (Fioc V Rocca-Volmerange 
1997), and galaxev (Bruzual V Chariot 2003). All these 
codes perform essentially the same computation. One adopts 
a star formation history (SFH) and an initial mass func¬ 
tion (IMF) to determine the present-day distribution of stel¬ 
lar masses and ages. Next, using a set of stellar evolution¬ 
ary tracks and atmospheres that give the luminosity (either 
frequency-dependent or integrated over some hlter) for a 
star of a particular mass and age, one integrates the stel¬ 
lar luminosity weighted by the age and mass distributions. 
These codes differ in the range of functional forms they al- 
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low for the IMF and SFH, and the evolutionary tracks and 
model atmospheres they use, but the underlying computa¬ 
tion is much the same in all of them. While this approach is 
adequate for many applications, it fails for systems with low 
masses and star formation rates (SFRs) because it implic¬ 
itly assumes that the stellar mass and age distributions are 
well-sampled. This is a very poor assumption both in star¬ 
forming dwarf galaxies and in resolved segments of larger 
galaxies. 

Significantly less work has been done in extending SPS 
methods to the regime where the IMF and SFH are not 
well-sampled. There are a number of codes available for 
simulating a simple stellar population (i.e., one where all 
the stars are the same age, so the SFH is described by a 
5 distribution) where the IMF is not well sampled (Mafz 
Apellaniz 2009; Popescu & Hanson 2009, 2010b, a; Foues- 
neau & Langon 2010; Fouesneau et al. 2012, 2014; Anders 
et al. 2013; de Meulenaer et al. 2013, 2014, 2015), and a great 
deal of analytic work has also been performed on this topic 
(Cerviho & Valls-Gabaud 2003; Cerviho & Luridiana 2004, 
2006) - see Cerviho (2013) for a recent review. However, 
these codes only address the problem of stochastic sampling 
of the IMF; for non-simple stellar populations, stochastic 
sampling of the SFH proves to be a more important effect 
(Fumagalli et al. 2011; da Silva et al. 2014). 

To handle this problem, we introduced the stochastic 
SPS code slug (da Silva et al. 2012), which includes full 
stochasticity in both the IMF and the SFH. Crucially, slug 
properly handles the clustered nature of star formation (e.g., 
Lada & Lada 2003; Krumholz 2014). This has two effects. 
First, clustering itself can interact with IMF sampling so 
that the total mass distribution produced by a clustered 
population is not identical to that of a non-clustered popu¬ 
lation drawn from the same underlying IMF. Second and 
perhaps more importantly, clustering causes large short- 
timescale fluctuations in the SFR even in galaxies whose 
mean SFR averaged over longer timescales is constant. Since 
its introduction, this code has been used in a number of ap¬ 
plications, including explaining the observed deficit of Ha 
emission relative to FUV emission in dwarf galaxies (Fu¬ 
magalli et al. 2011; Andrews et al. 2013, 2014), quantify¬ 
ing the stochastic uncertainties in SFR indicators (da Silva 
et al. 2014), analyzing the properties of star clusters in dwarf 
galaxies in the ANCST survey (Cook et al. 2012), and ana¬ 
lyzing Lyman continuum radiation from high-redshift dwarf 
galaxies (Forero-Romero & Dijkstra 2013). The need for a 
code with stochastic capabilities is likely to increase in the 
future, as studies of local galaxies such as PHAT (Dalcan- 
ton et al. 2012), HERACLES (Leroy et al. 2013), and LE- 
CUS (Calzetti et al. 2015), and even studies in the high red- 
shift universe (e.g., Jones et al. 2013a,b) increasingly push 
to smaller spatial scales and lower rates of star formation, 
where stochastic effects become increasingly important. 

In this paper we describe a major upgrade and expan¬ 
sion of slug, intended to make it a general-purpose solution 
for the analysis of stochastic stellar populations. This new 
version of slug allows essentially arbitrary functional forms 
for both the IME and the SEH, allows a wide range of stellar 
evolutionary tracks and atmosphere models, and can output 
both spectroscopy and photometry in a list of > 100 filters. 
It can also include the effects of reprocessing of the light 
by interstellar gas and by stochastically-varying amounts 


of dust, and can interface with the cloudy photoioniza¬ 
tion code (Eerland et al. 2013) to produce predictions for 
stochastically-varying nebular emission. Einally, we have 
coupled slug to a new set of tools for solving the “inverse 
problem” in stochastic stellar population synthesis: given a 
set of observed photometric properties, infer the posterior 
probability distribution for the properties of the underlying 
stellar population, in the case where the mapping between 
such properties and the photometry is stochastic and there¬ 
fore non-deterministic (e.g. da Silva et al. 2014). The full 
software suite is released under the GNU Public License, 
and is freely available from http: //www. slugsps. com.^ 

In the remainder of this paper, we describe slug and 
its companion software tools in detail (Section 2), and then 
provide a series of demonstrations of the capabilities of the 
upgraded version of the code (Section 3). We end with a 
summary and discussion of future prospects for this work 
(Section 4). 


2 THE SLUG SOFTWARE SUITE 

The slug software suite is a collection of tools designed to 
solve two related problems. The first is to determine the 
probability distribution function (PDF) of observable quan¬ 
tities (spectra, photometry, etc.) that are produced by a 
stellar population characterized by a specified set of phys¬ 
ical parameters (IMF, SFH, cluster mass function, and an 
array of additional ancillary inputs). This problem is ad¬ 
dressed by the core slug code (Section 2.1) and its adjunct 
cloudy_slug (Section 2.2), which perform Monte Carlo sim¬ 
ulations to calculate the distribution of observables. The 
second problem is to use those PDFs to solve the inverse 
problem: given a set of observed properties, what should we 
infer about the physical properties of the stellar population 
producing those observables? The bayesphot package pro¬ 
vides a general solution to this problem (Section 2.3), and 
the SFR_slug and cluster.slug packages specialize this gen¬ 
eral solution to the problem of inferring star formation rates 
for continuously star-forming galaxies, and masses, ages, 
and extinctions from simple stellar populations, respectively 
(Section 2.4). The entire software suite is available from 
http : //www. slugsps. com, and extensive documentation is 
available at http: //slug2 . readthedocs . org/en/latest/. 

2.1 slug: A Highly Flexible Tool for Stochastic 
Stellar Population Synthesis 

The core of the software suite is the stochastic stellar popula¬ 
tion synthesis code slug. The original slug code is described 
in da Silva et al. (2012), but the code described here is a com¬ 
plete re-implementation with greatly expanded capabilities. 
The code operates in two main steps; first it generates the 
spectra from the stellar population (Section 2.1.1) and then 
it post-processes the spectra to include the effects of nebular 
emission and dust absorption, and to produce photometric 

^ As of this writing http: //www. slugsps. com is hosted on Google 
Sites, and thus is inaccessible from mainland China. Chinese 
users can access the source code from https: //bitbucket. org/ 
krumholz/slug2, and should contact the authors by email for the 
ancillary data products. 
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predictions as well as spectra (Section 2.1.2). In this section 
we limit ourselves to a top-level description of the physi¬ 
cal model, and provide some more details on the numerical 
implementation in Appendices A, B, and C. 

2.1.1 Generating the Stellar Population 

Slug can simulate both simple stellar populations (SSPs, 
i.e., ones where all the stars are coeval) and composite ones. 
We begin by describing the SSP model, since composite stel¬ 
lar populations in slug are built from collections of SSPs. 
For an SSP, the stellar population is described by an age, a 
total mass, and an IMF. Slug allows nearly-arbitrary func¬ 
tional forms for the IMF, and is not restricted to a set of 
predetermined choices; see Appendix AI for details. 

In non-stochastic SPS codes, once an age, mass, and 
IMF are chosen, the standard procedure is to use a set 
of stellar tracks and atmospheres to predict the luminos¬ 
ity (either frequency-dependent, or integrated over a spec¬ 
ified filter) as a function of stellar mass and age, and to 
integrate the mass-dependent luminosity multiplied by the 
stellar mass distribution to compute the output spectrum 
at a specified age. Slug adds an extra step: instead of eval¬ 
uating an integral, it directly draws stars from the IMF un¬ 
til the desired total mass has been reached. As emphasized 
by a number of previous authors (e.g. Weidner & Kroupa 
2006; Haas & Anders 2010; Cerviho 2013; Popescu & Han¬ 
son 2014), when drawing a target mass Mtarget rather than 
a specified number of objects from a PDF, one must also 
choose a sampling method to handle the fact that, in gen¬ 
eral, it will not be possible to hit the target mass exactly. 
Many sampling procedures have been explored in the lit¬ 
erature, and slug provides a large number of options, as 
described in Appendix AI. 

Once a set of stars is chosen, slug proceeds much like 
a conventional SPS code. It uses a chosen set of tracks and 
atmospheres to determine the luminosity of each star, either 
wavelength-dependent or integrated over one or more pho¬ 
tometric filters, at the specified age. It then sums over all 
stars to determine the integrated light output. Details of the 
available sets of tracks and atmospheres, and slug’s method 
for interpolating them, are provided in Appendix A2. 

Composite stellar populations in slug consist of a col¬ 
lection of “star clusters”, each consisting of an SSP, plus 
a collection of “field stars” whose ages are distributed con¬ 
tinuously. In practice, this population is described by four 
parameters beyond those that describe SSPs: the fraction of 
stars formed in clusters as opposed to the field /c, the clus¬ 
ter mass function (CMF), cluster lifetime function (CLF), 
and star formation history (SFH). As with the IMF, the 
latter three are distributions, which can be specified using 
nearly arbitrary functional forms and a wide range of sam¬ 
pling methods as described in Appendix AI. For a simu¬ 
lation with a composite stellar population, slug performs 
a series of steps: (I) at each user-requested output time^, 

^ Slug can either output results at specified times, or the user 
can specify a distribution from which the output time is to be 
drawn. The latter capability is useful when we wish to sample a 
distribution of times continuously so that the resulting data set 
can be used to infer ages from observations - see Section 2.4. 


slug uses the SFH and the cluster fraction fc to compute 
the additional stellar mass expected to have formed in clus¬ 
ters and out of clusters (in “the field”) since the previous 
output; (2) for the non-clustered portion of the star forma¬ 
tion, slug draws the appropriate mass in individual stars 
from the IMF, while for the clustered portion it draws a 
cluster mass from the CMF, and then it fills each cluster 
with stars drawn from the IMF; (3) each star and star clus¬ 
ter formed is assigned an age between the current output 
and the previous one, selected to produce a realisation of 
the input SFH^; (4) each cluster that is formed is assigned 
a lifetime drawn from the CLF, which is the time at which 
the cluster is considered dispersed. 

The end result of this procedure is that, at each output 
time, slug has constructed a “galaxy” consisting of a set of 
star clusters and field stars, each with a known age. Once 
the stellar population has been assembled, the procedure 
for determining spectra and photometry is simply a sum 
over the various individual simple populations. In addition 
to computing the integrated spectrum of the entire popula¬ 
tion, slug can also report individual spectra for each cluster 
that has not disrupted (i.e., where the cluster age is less 
than the value drawn from the CLF for that cluster). Thus 
the primary output consists of an integrated monochromatic 
luminosity Lx for the entire stellar population, and a value 
Lx,i for the zth remaining cluster, at each output time. 

Finally, we note that slug is also capable of skipping 
the sampling procedure and evaluating the light output of 
stellar populations by integrating over the IMF and SFH, 
exactly as in a conventional, non-stochastic SPS code. In 
this case, slug essentially emulates starburst99 (Leitherer 
et al. 1999; Vazquez V Leitherer 2005), except for sub¬ 
tle differences in the interpolation and numerical integra¬ 
tion schemes used (see Appendix A2). Slug can also behave 
“semi-stochastically”, evaluating stars’ contribution to the 
light output using integrals to handle the IMF up to some 
mass, and using stochastic sampling at higher masses. 


2.1.2 Post-Proeessing the Speetra 

Once slug has computed a spectrum for a simple or com¬ 
posite stellar population, it can perform three additional 
post-processing steps. First, it can provide an approximate 
spectrum after the starlight has been processed by the sur¬ 
rounding H II region. In this case, the nebular spectrum is 
computed for an isothermal, constant-density H ll region, 
within which it is assumed that He is all singly-ionized. Un¬ 
der these circumstances, the photoionized volume U, elec¬ 
tron density Ue, and hydrogen density nn obey the relation 


^ One subtlety to note here is that the choice of output grid can 
produce non-trivial modifications of the statistical properties of 
the output in cases where the expected mass of stars formed is 
much smaller than the maximum possible cluster mass. For exam¬ 
ple, consider a simulation with a constant SFR and an output grid 
chosen such that the expected mass of stars formed per output 
time is 10^ Mq. If the sampling method chosen is ST0P_NEAREST 
(see Appendix AI), the number of 10® Mq clusters formed will 
typically be smaller than if the same SFH were used but the out¬ 
puts were spaced 100 times further apart, giving an expected mass 
of 10® Mq between outputs. 
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is the hydrogen-ionizing luminosity, Ljy = X^Lxjc is the lu¬ 
minosity per unit frequency, /(H°) = 13.6 eV is the ioniza¬ 
tion potential of neutral hydrogen, Ue = (1 + XHe)^H, xne is 
the He abundance relative to hydrogen (assumed to be 0.1), 
(j) is the fraction of H-ionizing photons that are absorbed by 
hydrogen atoms within the H ii region rather than escaping 
or being absorbed by dust, and a^^\T) is the temperature- 
dependent case B recombination coefficient. The nebular lu¬ 
minosity can then be written as 


TA,neb — 'Tneb^e^nH — ^neb^ 


Q(H°) 

a(B)(T)’ 


(3) 


where 7neb is the wavelength-dependent nebular emission 
coefficient. Our calculation of this quantity includes the fol¬ 
lowing processes: free-free and bound-free emission arising 
from electrons interacting with and He^, two-photon 
emission from neutral hydrogen in the 2s state, H recombi¬ 
nation lines, and non-H lines computed approximately from 
tabulated values. Full details on the method by which we 
perform this calculation are given in Appendix B. The com¬ 
posite spectrum after nebular processing is zero at wave¬ 
length shorter than 912 A (corresponding to the ionization 
potential of hydrogen), and the sum of intrinsic stellar spec¬ 
trum L\ and the nebular spectrum La, neb at longer wave¬ 
lengths. Slug reports this quantity both cluster-by-cluster 
and for the galaxy as a whole. 

Note that the treatment of nebular emission included in 
slug is optimized for speed rather than accuracy, and will 
produce signihcantly less accurate results than cloudy_slug 
(see Section 2.2). The major limitations of the built-in ap¬ 
proach are: (1) because slug’s built-in calculation relies on 
pre-tabulated values for the metal line emission and as¬ 
sumes either a constant or a similarly pre-tabulated temper¬ 
ature (which signihcantly affects the continuum emission), it 
misses the variation in emission caused by the fact that H ii 
regions exist at a range of densities and ionization param¬ 
eters, which in turn induces substantial variations in their 
nebular emission (e.g. Yeh et al. 2013; Verdolini et al. 2013); 
(2) slug’s built-in calculation assumes a uniform density H ii 
region, an assumption that will fail at high ionizing lumi¬ 
nosities and densities due to the effects of radiation pressure 
(Dopita et al. 2002; Krumholz et al. 2009; Draine 2011; Yeh 
& Matzner 2012; Yeh et al. 2013); (3) slug’s calculation 
correctly captures the effects of stochastic variation in the 
total ionizing luminosity, but it does not capture the effects 
of variation in the spectral shape of the ionizing continuum, 
which preliminary testing suggests can cause variations in 
line luminosities at the few tenths of a dex level even for 
hxed total ionizing hux. The reason for accepting these lim¬ 
itations is that the assumptions that cause them also make 
it possible to express the nebular emission in terms of a few 
simple parameter that can be pre-tabulated, reducing the 
computational cost of evaluating the nebular emission by 
orders of magnitude compared to a fully accurate calcula¬ 
tion with cloudy_slug. 

The second post-processing step available is that slug 
can apply extinction in order to report an extincted spec¬ 


trum, both for the pure stellar spectrum and for the nebula- 
processed spectrum. Extinction can be either hxed to a con¬ 
stant value for an entire simulation, or it can be drawn from 
a specihed PDF. In the latter case, for simulations of com¬ 
posite stellar populations, every cluster will have a different 
extinction. More details are given in Appendix B. 

As a third and hnal post-processing step, slug can con¬ 
volve all the spectra it computes with one or more hlter 
response functions in order to predict photometry. Slug in¬ 
cludes the large list of hlter response functions maintained 
as part of FSPS (Conroy et al. 2010; Conroy & Gunn 2010), 
as well as a number of Hubble Space Telescope hlters^ not 
included in FSPS; at present, more than 130 hlters are avail¬ 
able. As part of this calculation, slug can also output the 
bolometric luminosity, and the photon luminosity in the H°, 
He°, and He+ ionizing continua. 


2.2 cloudy_slug: Stochastic Nebular Line Emission 

In broad outlines, the cloudy_slug package is a software in¬ 
terface that automatically takes spectra generated by slug 
and uses them as inputs to the cloudy radiative transfer 
code (Ferland et al. 2013). Cloudy_slug then extracts the 
continuous spectra and lines returned by cloudy, convolves 
them with the same set of hlters as in the original slug 
calculation in order to predict photometry. The user can 
optionally also examine all of cloudy’s detailed outputs. As 
with the core slug code, details on the software implemen¬ 
tation are given in Appendix C. 

When dealing with composite stellar populations, cal¬ 
culating the nebular emission produced by a stellar pop¬ 
ulation requires making some physical assumptions about 
how the stars are arranged, and cloudy_slug allows two ex¬ 
treme choices that should bracket reality. One extreme is to 
assume that all stars are concentrated in a single point of 
ionizing radiation at the center of a single H ii region. We 
refer to this as integrated mode, and in this mode the only 
free parameters to be specihed by a user are the chemical 
composition of the gas into which the radiation propagates 
and its starting density. 

The opposite assumption, which we refer to as cluster 
mode, is that every star cluster is surrounded by its own 
H II region, and that the properties of these regions are 
to be computed individually using the stellar spectrum of 
each driving cluster and only then summed to produce a 
composite output spectrum. At present cluster mode does 
not consider contributions to the nebular emission from held 
stars, and so should not be used when the cluster fraction 
fc < 1. In the cluster case, the H il regions need not all 
have the same density or radius; indeed, one would expect 
a range of both properties to be present, since not all star 
clusters have the same age or ionizing luminosity. We handle 
this case using a simplihed version of the H ll population 
synthesis model introduced by Verdolini et al. (2013) and 
Yeh et al. (2013). The free parameters to be specihed in this 
model are the chemical composition and the density in the 


^ These hlters were downloaded from http://www. 
stsci.edu/~WFC3/UVIS/SystemThroughput/ and http: 
//www.stsci.edu/hst/acs/analysis/throughputs for the 
UVIS and ACS instruments, respectively. 
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ambient neutral ISM around each cluster. For an initially- 
neutral region of uniform hydrogen number density nu , the 
radius of the H ll at a time t after it begins expanding is 
well-approximated by (Krumholz & Matzner 2009) 
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where = 2.59 x 10“^^ cm^ s“^ is the case B recom¬ 
bination coefficient at 10^ K, Tn = 10^ is the typical H ii 
region temperature, /trap = 2 is a trapping factor that ac¬ 
counts for stellar wind and trapped infrared radiation pres¬ 
sure, 'ip = 3.2 is the mean photon energy in Rydberg for a 
fully sampled IMF at zero age^, and /i = 1.33 is the mean 
mass per hydrogen nucleus for gas with the usual helium 
abundance xue = 0.1. The solution includes the effects of 
both radiation and gas pressure in driving the expansion. 
Once the radius is known, the density near the H ll region 
center (the parameter required by cloudy) can be computed 
from the usual ionization balance argument. 


nil 


n \ 1/2 

3Q(H°) \ 
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Note that the factor of 4.4 in the denominator accounts 
for the extra free electrons provided by helium, assum¬ 
ing it is singly ionized. Also note that this will not give 
the correct density during the brief radiation pressure- 
dominated phase early on in the evolution, but that this 
period is very brief (though the effects of the extra boost 
provided by radiation pressure can last much longer), and 
no simple analytic approximation for this density is avail¬ 
able (Krumholz & Matzner 2009). Also note that, although 
cloudy_slug presently only supports this parameterization 
of H II region radius and density versus time, the code is sim¬ 
ply a python script. It would therefore be extremely straight¬ 
forward for a user who prefers a different parameterization 
to alter the script to supply it. 

In cluster mode, cloudy_slug uses the approximation 
described above to compute the density of the ionized gas 
in the vicinity of each star cluster, which is then passed as 
an input to cloudy along with the star cluster’s spectrum. 
Note that computations in cluster mode are much more ex¬ 
pensive than those in integrated mode, since the latter re¬ 
quires only a single call to cloudy per time step, while the 
former requires one per cluster. To ease the computational 


^ This value of is taken from Murray & Rahman (2010), and is 
based on a Chabrier (2005) IMF and a compilation of empirically- 
measured ionizing luminosities for young stars. However, alterna¬ 
tive IMFs and methods of computing the ionizing luminosity give 
results that agree to tens of percent. 


requirements slightly, in cluster mode one can set a thresh¬ 
old ionizing luminosity below which the contribution to the 
total nebular spectrum is assumed to be negligible, and is 
not computed. 


2.3 bayesphot: Bayesian Inference from Stellar 
Photometry 

2.3.1 Description of the Method 


Simulating stochastic stellar populations is useful, but to 
fully exploit this capability we must tackle the inverse prob¬ 
lem: given an observed set of stellar photometry, what 
should we infer about the physical properties of the un¬ 
derlying stellar population in the regime where the map¬ 
ping between physical and photometric properties is non- 
deterministic? A number of authors have presented numer¬ 
ical methods to tackle problems of this sort, mostly in 
the context of determining the properties of star clusters 
with stochastically-sampled IMFs (Popescu & Hanson 2009, 
2010b, a; Fouesneau & Langon 2010; Fouesneau et al. 2012; 
Popescu et al. 2012; Asa’d & Hanson 2012; Anders et al. 
2013; de Meulenaer et al. 2013, 2014, 2015), and in some 
cases in the context of determining SFRs from photometry 
(da Silva et al. 2014). Our method here is a generalization of 
that developed in da Silva et al. (2014), and has a number of 
advantages as compared to earlier methods, both in terms 
of its computational practicality and its generality. 

Consider stellar systems characterized by a set of N 
physical parameters x = (xi, 0 : 2 , .. .xn)] in the example of 
star clusters below we will have iV = 3, with xi, X 2 , and X 3 
representing the logarithm of the mass, the logarithm of the 
age, and the extinction, while for galaxies forming stars at 
a continuous rate we will have A/ = 1, with xi representing 
the logarithm of the SFR. The light output of these systems 
is known in M photometric bands; let y = (yi, 2 / 2 , • • • Vm) be 
a set of photometric values, for example magnitudes in some 
set of filters. Suppose that we observe a stellar population in 
these bands and measure a set of photometric values yobs, 
with some errors cry = {ay ^, ay 2 , • • •, which for sim¬ 

plicity we will assume are Gaussian. We wish to infer the 
posterior probability distribution for the physical parame¬ 
ters given the observed photometry and photometric errors, 
p(x I yobs;0-y)- 

Following da Silva et al. (2014), we compute the pos¬ 
terior probability via implied conditional regression coupled 
to kernel density estimation. Let the joint PDF of physical 
and photometric values for the population of all the stellar 
populations under consideration be p(x,y). We can write 
the posterior probability distribution we seek as 


p(x I y) oc 


P(x,y) 

p(y) 


( 11 ) 


where p{y) is the distribution of the photometric variables 
alone, i.e., 

P{y) J p{x,y)dx. (12) 

If we have an exact set of photometric measurements yobs, 
with no errors, then the denominator in equation (11) is 
simply a constant that will normalize out, and the posterior 
probability distribution we seek is distributed simply as 
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p(x I yobs) OC p(x, yobs). (13) 

In this case, the problem of computing p(x | yobs) reduces 
to that of computing the joint physical-photometric PDF 
p(x, y) at any given set of observed photometric values yobs- 
For the more general case where we do have errors, we 
first note that the posterior probability distribution for the 
true photometric value y is given by the prior probability 
distribution for photometric values multiplied by the likeli¬ 
hood function associated with our measurements. The prior 
PDF of photometric values is simply p(y) as given by equa¬ 
tion (12), and for a central observed value of yobs with errors 
(Ty, the likelihood function is simply a Gaussian. Thus the 
PDF of y given our observations is 


P(y I yobs) OC p(y)G(y - yobs, (Ty) 


(14) 


where 


G{y,cry) OC exp 


(y\ 


, v\x y 


+ +•■ 

2<^?m A 


(15) 


is the usual multi-dimensional Gaussian function. The pos¬ 
terior probability for the physical parameters is then simply 
the convolution of equations (11) and (14), i.e.. 


p(x I yobs) 


J p(x I y)p(y I yobs)dy 
J P(x,y) G(y -yobs,o-y)dy. 


(16) 


Note that we recover the case without errors, equation (13), 
in the limit where cry ^ 0, because in that case the Gaussian 
G(y - yobs, (Ty) (5(y - yobs)- 

We have therefore reduced the problem of computing 
p(x I yobs) to that of computing p(x, y), the joint PDF of 
physical and photometric parameters, for our set of stellar 
populations. To perform this calculation, we use slug to 
create a large library of models for the type of stellar pop¬ 
ulation in question. We draw the physical properties of the 
stellar populations (e.g., star cluster mass, age, and extinc¬ 
tion) in our library from a distribution pub (x), and for each 
stellar population in the library we have a set of physical 
and photometric parameters x^ and y^. We estimate p(x, y) 
using a kernel density estimation technique. Specifically, we 
approximate this PDF as 


P(x,y) OC y^Wiii:(x-Xi,y-yi;h), 
i 


(17) 


where wi is the weight we assign to sample z, iF(x; h) is our 
kernel function, h = (h^i ,hx 2 , - • •, , hy^, hy^,..., hy^) 

is our bandwidth parameter (see below), and the sum runs 
over every simulation in our library. We assign weights to en¬ 
sure that the distribution of physical parameters x matches 
whatever prior probability distributions we wish to assign 
for them. If we let Pprior(x) represent our priors, then this is 


Pprior (^i) 
Plib(Xi) 


(18) 


Note that, although we are free to choose piib(x) = Pprior(x) 
and thus weight all samples equally, it is often advantageous 
for numerical or astrophysical reasons not to do so, because 
then we can distribute our sample points in a way that is 
chosen to maximize our knowledge of the shape of p(x, y) 
with the fewest possible realizations. As a practical example, 
we know that photometric outputs will fluctuate more for 


galaxies with low star formation rates than for high ones, so 
the library we use for SFR_slug (see below) is weighted to 
contain more realizations at low than high SFR. 

We choose to use a Gaussian kernel function, i^(x; h) = 
G(x, h), because this presents significant computational ad¬ 
vantages. Specifically, with this choice, equation (16) be¬ 
comes 

p(x I yobs) 

OC y^ WiG{x - Xi,ha;) 

i 

j G(y - yobs, (Ty) G(y - y*, h^) dy 

OC y^ WiG{x - Xi,ha;) 

i 

G(yob8 - yi, \j(Tl + hg) 

= y^ WiGjjx - Xi,yob3 - yb,h') 

i 

where hx = {hx^, hx 2 , ■ ■ ■, hxj^) is the bandwidth in the 
physical dimensions, = (^yi? hy^, - • •, hyj^) is the band¬ 
width in the photometric dimensions, and the quadrature 
sum VH + is understood to be computed independently 
over every index in cry and h^. The new quantity we have 
introduced, 

h' = (hj;, +0-2), (22) 

is simply a modified bandwidth in which the bandwidth in 
the photometric dimensions has been broadened by adding 
the photometric errors in quadrature sum with the band- 
widths in the corresponding dimensions. Note that, in going 
from equation (19) to equation (20) we have invoked the re¬ 
sult that the convolution of two Gaussians is another Gaus¬ 
sian, whose width is the quadrature sum of the widths of 
the two input Gaussians, and whose center is located at the 
difference between the two input centers. 

As an adjunct to this result, we can also immediately 
write down the marginal probabilities for each of the phys¬ 
ical parameters in precisely the same form. The marginal 
posterior probability distribution for xi is simply 

p{xi I yobs) OC J p(x I yobs) dx 2 dx3 - - - dxN (23) 

OC y^^WiG{xi - xi^i, hi) 

i 

G(yobs-yi,A/o-^ + h2), (24) 

and similarly for all other physical variables. This expression 
also immediately generalizes to the case of joint marginal 
probabilities of the physical variables. We have therefore suc¬ 
ceeded in writing down the posterior probability distribution 
for all the physical properties, and the marginal posterior 
distributions of any of them individually or in combination, 
via a kernel density estimation identical to that given by 
equation (17). One advantage of equations (21) and (24) is 
that they are easy to evaluate quickly using standard nu¬ 
merical methods. Details on our software implementation 
are given in Appendix G. 


(19) 

( 20 ) 
( 21 ) 
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2.3.2 Error Estimation and Properties of the Library 


An important question for bayesphot, or any similar 
method, is how large a library must be in order to yield 
reliable results. Of course this depends on what quantity is 
being estimated and on the shape of the underlying distri¬ 
bution. Full, rigorous, error analysis is best accomplished by 
bootstrap resampling. However, we can, without resorting to 
such a computationally-intensive procedure, provide a use¬ 
ful rule of thumb for how large a library must be so that 
the measurements rather than the library are the dominant 
source of uncertainty in the resulting derived properties. 

Consider a library of n simulations, from which we wish 
to measure the value Xq that delineates a certain quantile 
(i.e., the percentile p = lOOg) of the underlying distribution 
from which the samples are drawn. Let Xi for z = 1... n be 
the samples in our library, ordered from smallest to largest. 
Our central estimate for the value of Xq that delineates quin¬ 
tile q, will naturally be Xqn- (Throughout this argument we 
will assume that n is large enough that we need not worry 
about the fact that ranks in the list, such as gn, will not 
exactly be integers. Extending this argument to include the 
necessary continuity correction adds mathematical compli¬ 
cation but does not change the basic result.) For example, 
we may have a library of 10® simulations of galaxies at a par¬ 
ticular SFR (or in a particular interval of SFRs), and wish 
to know what ionizing luminosity corresponds to the 95th 
percentile at that SFR. Our estimate for the 95th percentile 
ionizing luminosity will simply be the ionizing luminosity of 
the 950,000th simulation in the library. 

To place confidence intervals around this estimate, note 
that the probability that a randomly chosen member of the 
population will have a value x < Xq is q, and conversely the 
probability that x > Xq is 1 — q. Thus in our library of n 
simulations, the probability that we have exactly k samples 
for which x < Xq is given by the binomial distribution. 


When n ^ 1, k ^ 1, and n — k ^ 1, the binomial distribu¬ 
tion approaches a Gaussian distribution of central value qn 
and standard deviation q(l — q)n: 


Pr(/c) 


1 

i/27rg(l - q)n 


exp 


{k — qnY 

2q(l - q)n 


(26) 


That is, the position n in our library of simulations such that 
Xi < Xq for alH < n, and Xi > Xq for all z > n, is distributed 
approximately as a Gaussian, centered at qn, with disper¬ 
sion q(l — q)n. This makes it easy to compute confidence 
intervals on Xq, because it reduces the problem that of com¬ 
puting confidence intervals on a Gaussian. Specifically, if we 
wish to compute the central value and confidence interval c 
(e.g., c = 0.68 is the 68% confidence interval) for the rank r 
corresponding to percentile q, this is 

r ~ ± y/2^(r^^^)nerf~^(c). (27) 


The corresponding central value for Xq (as opposed to its 
rank in the list) is Xgn, and the confidence interval is 

^\f^qn — ^JEq{E^q)ner^~^ (c) ’ ^qn-\-^f 2q{l — q)n erf ~ ^ (c) j ^ 

In our example of determining the 95th percentile from a 
library of 10® simulations, our central estimate of this value 


is the X 95 o,ooo, the 950,000th simulation in the library, and 
the 90% confidence interval is 359 ranks on either side of this. 
That is, our 90% confidence interval extends from 2 : 949,641 
to X950,359- 

This method may be used to define confidence intervals 
on any quantile derived from a library of slug simulations. 
Of course this does not specifically give confidence intervals 
on the results derived when such a library is used as the 
basis of a Bayesian inference from bayesphot. However, this 
analysis can still provide a useful constraint: the error on 
the values derived from bayesphot cannot be less than the 
errors on the underlying photometric distribution we have 
derived from the library. Thus if we have a photometric mea¬ 
surement that lies at the qih quantile of the library we are 
using for kernel density estimation in bayesphot, the ef¬ 
fective uncertainty our kernel density estimation procedure 
provides is at a minimum given by the uncertainty in the 
quantile value Xq calculated via the procedure above. If this 
uncertainty is larger than the photometric uncertainties in 
the measurement, then the resolution of the library rather 
than the accuracy of the measurement will be the dominant 
source of error. 


2.4 sfr_slug and cluster_slug: Bayesian Star 

Formation Rates and Star Cluster Properties 

The slug package ships with two Bayesian inference mod¬ 
ules based on bayesphot. SFR_slug, first described in da 
Silva et al. (2014) and modified slightly to use the improved 
computational technique described in the previous section, 
is a module that infers the posterior probability distribution 
of the SFR given an observed flux in Ho, GALEX FUV, or 
bolometric luminosity, or any combination of the three. The 
library of models on which it operates (which is included in 
the slug package) consists of approximately 1.8 x 10® galax¬ 
ies with constant SFRs sampling a range from 10~® — 10® ® 
Mq yr“^, with no extinction; the default bandwidth is 0.1 
dex. We refer readers for da Silva et al. (2014) for a further 
description of the model library. SFR_slug can download the 
default library automatically, and it is also available as a 
standalone download from http : //www. slugsps . com/data. 

The cluster_slug package performs a similar task for 
inferring the mass, age, and extinction of star clusters. The 
slug models that cluster_slug uses consist of libraries of 
10^ simple stellar populations with masses in the range 
log(M/Mo) = 2-8, ages log(T/yr) = 5 - logTmax, and 
extinctions Ay = 0 — 3 mag. The maximum age Tmax is 
either 1 Gyr or 15 Gyr, depending on the choice of tracks 
(see below). The data are sampled uniformly in Ay - In mass 
the library sampling is dN/dM oc for masses up to 10^ 
Mq, and as dN/dM oc M~^ at higher masses; similarly, the 
library age sampling is dN/dT oc T~^ for T < 10® yr, and as 
dN/dT oc T~^ at older ages. The motivation for this sam¬ 
pling is that it puts more computational effort at younger 
ages and low masses, where stochastic variation is greatest. 
The libraries all use a Ghabrier (2005) IMF, and include 
nebular emission computed using 0 = 0.73 and an ionization 
parameter log lA — —3 (see Appendix B). Libraries are avail¬ 
able using either Padova or Geneva tracks, with the former 
having a maximum age of 15 Gyr and the latter a maximum 
age of 1 Gyr. The Geneva tracks are available using either 
Milky Way or “starburst” extinction curves (see Appendix 
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Stop After 
Stop Nearest 
Stop Before 
Sorted Sampling 



Figure 1. Histograms of the mass of the most massive stars in 
sets of 10"^ “cluster” simulations, which were performed with dif¬ 
ferent sampling techniques (top panel) and different target cluster 
masses (bottom panel). In the top panel, four different stop crite¬ 
ria are adopted to sample a Kroupa (2002) IMF for clusters of 50 
Mq. In the bottom panel, the default ST0P_NEAREST condition is 
chosen to sample from a Kroupa (2002) IMF in clusters of three 
different masses. 


B). The default bandwidth is 0.1 dex in mass and age, 0.1 
mag in extinction, and 0.25 mag (corresponding to 0.1 dex 
in luminosity) in photometry. For each model, photometry 
is computed for a large range of filters, listed in Table 1. 
As with SFR_slug, cluster_slug can automatically down¬ 
load the default library, and the data are also available as a 
standalone download from http: //www. slugsps . com/data. 
Full spectra for the library, allowing the addition of further 
filters as needed, are also available upon request; they are 
not provided for web download due to the large file sizes 
involved. 


3 SAMPLE APPLICATIONS 

In this section, we present a suite of test problems with the 
goal of illustrating the different capabilities of slug, high¬ 
lighting the effects of key parameters on slug’s simulations, 
and validating the code. Unless otherwise stated, all com¬ 
putations use the Geneva (2013) non-rotating stellar tracks 
and stellar atmospheres following the star burst 99 imple¬ 
mentation. 


3.1 Sampling Techniques 

As discussed above and in the literature (e.g., Weidner & 
Kroupa 2006; Haas & Anders 2010; da Silva et al. 2012; 
Cerviho et al. 2013; Popescu & Hanson 2014), the choice 
of sampling technique used to generate stellar masses can 
have signihcant effects on the distribution of stellar masses 


and the hnal light output, even when the underlying distri¬ 
bution being sampled is held hxed. This can have profound 
astrophysical implications. Variations in sampling technique 
even for a hxed IMF can produce systematic variations in 
galaxy colors (e.g., Haas V Anders 2010), nucleosynthetic 
element yields (e.g., Koppen et al. 2007; Haas V Anders 
2010 ), and observational tracers of the star formation rate 
(e.g., Weidner et al. 2004; Phamm-Altenburg et al. 2007). It 
is therefore of interest to explore how sampling techniques 
influence various aspects of stellar populations. We do so as 
a hrst demonstration of slug’s capabilities. Here we con¬ 
sider the problem of populating clusters with stars from an 
IMF, but our discussion also applies to the case of “galaxy” 
simulations. Specihcally, we run four sets of 10^ “cluster” 
simulations with a target mass of 50 M© by sampling a 
Kroupa (2002) IMF with the ST0P_NEAREST, ST0P_AFTER, 
ST0P_BEF0RE, SORTED.SAMPLING conditions (see Appendix 
Al). In the following, we analyse a single timestep at 10® 
yr. 

By default, slug adopts the ST0P_NEAREST condition, 
according to which the hnal draw from the IMF is added to 
the cluster only if the inclusion of this last star minimises the 
absolute error between the target and the achieved cluster 
mass. In this case, the cluster mass sometimes exceeds and 
sometimes falls short of the desired mass. The STOP .AFTER 
condition, instead, always includes the hnal draw from the 
IMF. Thus, with this choice, slug simulations produce clus¬ 
ters with masses that are always in excess of the target 
mass. The opposite behaviour is obviously recovered by the 
ST0P.BEF0RE condition, in which the hnal draw is always 
rejected. Finally, for SORTED .SAMPLING condition, the hnal 
cluster mass depends on the details of the chosen IMF. 

Besides this manifest effect of the sampling techniques 
on the achieved cluster masses, the choice of sampling pro¬ 
duces a drastic effect on the distribution of stellar masses, 
even for a hxed IMF. This is shown in the top panel of 
Figure 1, where we display histograms for the mass of the 
most massive stars within these simulated clusters. One 
can see that, compared to the default STOP .NEAREST con¬ 
dition, the ST0P_AFTER condition results in more massive 
stars being included in the simulated clusters. Conversely, 
the ST0P_BEF0RE and the SORTED .SAMPLING undersample the 
massive end of the IMF. Such a different stellar mass distri¬ 
bution has direct implications for the photometric proper¬ 
ties of the simulated stellar populations, especially for wave¬ 
lengths that are sensitive to the presence of most massive 
stars. Similar results have previously been obtained by other 
authors, including Haas V Anders (2010) and Cerviho et al. 
(2013), but prior to slug no publicly-available code has been 
able to tackle this problem. 

The observed behaviour on the stellar mass distribu¬ 
tion for a particular choice of sampling stems from a generic 
mass constraint: clusters cannot be hhed with stars that are 
more massive than the target cluster itself. Slug does not en¬ 
force this condition strictly, allowing for realizations in which 
stars more massive than the target cluster mass are included. 
However, different choices of the sampling technique result 
in a different degree with which slug enforces this mass con¬ 
straint. To further illustrate the relevance of this effect, we 
run three additional “cluster” simulations assuming again a 
Kroupa (2002) IMF and the default ST0P.NEAREST condition. 
Each simulation is composed by 10^ trials, but with a target 
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Table 1. Filters in the cluster .slug library 


Type 

Filter Name 

HST WFC3 UVIS wide 

F225W, F275W, F336W, F360W, F438W, F475W, F555W, F606W, F775W, F814W 

HST WFC3 UVIS medium/narrow 

F547M, F657N, F658N 

HST WFC3 IR 

F098M, F105W, FllOW, F125W, F140W, F160W 

HST ACS wide 

F435W, F475W, F555W, F606W, F625W, F775W, F814W 

HST ACS medium/narrow 

F550M, F658N, F660N 

HST ACS HRC 

F330W 

HST ACS SBC 

F125LP, F140LP, F150LP F165LP 

Spitzer IRAC 

3.6, 4.5, 5.8, 8.0 

GALEX 

FUV, NUV 

Johnson-Cousins 

U, B, V, R, I 

SDSS 

u, g, r, i, z 

2MASS 

J, H, Ks 


cluster mass of Mci,t = 50, 500, and 10^ M©. The bottom 
panel of Figure 1 shows again the mass distribution of the 
most massive star in each cluster. As expected, while massive 
stars cannot be found in low mass clusters, the probability 
of finding at least a star as massive as the IMF upper limit 
increases with the cluster mass. At the limit of very large 
masses, a nearly fully-sampled IMF is recovered as the mass 
constraint becomes almost irrelevant. Again, similar results 
have previously been obtained by Cervino (2013). 

3.2 Stochastic Spectra and Photometry for 
Varying IMFs 

Together with the adopted sampling technique, the choice 
of the IMF is an important parameter that shapes slug sim¬ 
ulations. We therefore continue the demonstration of slug’s 
capabilities by computing spectra and photometry for simple 
stellar populations with total masses of 500 M© for three dif¬ 
ferent IMFs. We create 1000 realizations each of such a pop¬ 
ulation at times from 1 — 10 Myr in intervals of 1 Myr, using 
the IMFs of Chabrier (2005), Kroupa (2002), and Weidner 
Sz Kroupa (2006). The former two IMFs use ST0P_NEAREST 
sampling, while the latter uses S0RTED_SAMPLING, and is oth¬ 
erwise identical to the Kroupa (2002) IMF. 

Figures 2-4 show the distributions of spectra and 
photometry that result from these simulations. Stochastic 
variations in photometry have been studied previously by a 
number of authors, going back to Chiosi et al. (1988), but to 
our knowledge no previous authors have investigated simi¬ 
lar variations in spectra. The plots also demonstrate slug’s 
ability to evaluate the full probability distribution for both 
spectra and photometric filters, and reveal interesting phe¬ 
nomena that would not be accessible to a non-stochastic SPS 
code. In particular, Figure 2 shows that the mean specific 
luminosity can be orders of magnitude larger than the mode. 
The divergence is greatest at wavelengths of a few hundred 
A at ages ~ 2 — 4 Myr, and wavelengths longer than ~ 5000 
A at 10 Myr. Indeed, at 4 Myr, it is noteworthy that the 
mean spectrum is actually outside the 10 — 90th percentile 
range. In this particular example, the range of wavelengths 
at 4 Myr where the mean spectrum is outside the 10 — 90th 
percentile corresponds to energies of ~ 3 Ryd. For a stellar 
population 4 Myr old, these photons are produced only by 
WR stars, and most prolifically by extremely hot WC stars. 
It turns out that a WC star is present ~ 5% of the time. 



A [A] 


Figure 3. Photometry of simple stellar populations at the same 
times as shown in Figure 2. In each panel, the thin lines show the 
same mean spectra plotted in Figure 2. Circles with error bars 
show the specific luminosity Lx measured in each of the indicated 
filters: GALEX FUV, and HST UVIS F225W, F336W, F555W, 
and F814W. The abcissae of the green points are placed at the 
effective wavelength of each filter, with the red and blue offset to 
either side for clarity. The labeled black curves in the top panel 
show the filter response functions for each filter. For these points, 
filled circles indicate the mean value, open circles indicate the 
median value, and error bars indicate the range from the 10 —90th 
percentile. The leftmost, square points with error bars show the 
comparable mean, median, and range for the ionizing luminosity 
(scale on the right axis). 

but these cases are so much more luminous than when a 
WC star is not present that they are sufficient to drag the 
mean upward above the highest luminosities that occur in 
the ~ 95% of cases when no WC star is present. 

A similar phenomenon is visible in the photometry 
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Figure 2. Spectra of simple stellar populations at ages of 1,2,4 and 10 Myr (top to bottom rows) for the IMFs of Chabrier (2005, 
left column), Kroupa (2002, middle column), and Weidner & Kroupa (2006, left column). In each panel, thick black lines indicate the 
mean, and points show the locations of individual simulations (10% of the simulations, selected at random), with the color of the point 
indicating the probability density for the monochromatic luminosity at each wavelength (see color bar) Probability densities are evaluated 
independently at each monochromatic luminosity, and are normalized to have a maximum of unity at each wavelength.. 


shown by Figure 3. In most of the filters the 10 — 90th per¬ 
centile range is an order of magnitude wide, and for the 
ionizing luminosity and the HST UVIS F814W at 10 Myr 
the spread is more than two orders of magnitude, with signif¬ 
icant offsets between mean and median indicating a highly 
asymmetric distribution. Figure 4, which shows the full dis¬ 
tributions for several of the filters, confirms this impres¬ 
sion: at early times the cumulative distribution functions 
are extremely broad, and the ionizing luminosity in partic¬ 
ular shows a broad distribution at low Q(H°) and then a 
small number of simulations with large (5(H°). Note that 
the percentile values are generally very well-determined by 
our set of 1000 simulations. Using the method described in 
Section 2.3.2 to compute the errors on the percentiles, we 
find that the 68% confidence interval on the 10th, 50th, and 
90th percentile values is less than 0.1 dex wide at almost 
all times, filters, and IMFs. The only exception is at 1 Myr, 
where the 68% confidence interval on the 10th percentile of 
ionizing luminosity is ~ 0.2 — 0.3 dex wide. 

The figures also illustrate slug’s ability to capture the 
“IGIMF effect” (Weidner & Kroupa 2006) whereby sorted 
sampling produces a systematic bias toward lower luminosi¬ 
ties at short wavelengths and early times. Both the spectra 
and photometric values for the Weidner & Kroupa (2006) 
IMF are systematically suppressed relative to the IMFs that 
use a sampling method that is less biased against high mass 
stars (cf. Figure 1). 


3.3 Cluster Fraction and Cluster Mass Function 

The first two examples have focused on simple stellar popu¬ 
lations, namely collections of stars that are formed at coeval 
times to compose a “cluster” simulation. Our third example 
highlights slug’s ability to simulate composite stellar pop¬ 
ulations. Due to slug’s treatment of stellar clustering and 
stochasticity, simulations of such populations in slug differ 
substantially from those in non-stochastic SPS codes. Unlike 
in a conventional SPS code, the outcome of a slug simula¬ 
tion is sensitive to both the CMF and the fraction fc of stars 
formed in clusters, slug can therefore be applied to a wide 
variety of problems in which the fraction of stars formed 
within clusters or in the field becomes a critical parameter 
(e.g. Fumagalli et al. 2011). 

The relevance of clusters in slug simulations arises from 
two channels. First, because stars form in clusters of finite 
size, and the interval between cluster formation events is 
not necessarily large compared to the lifetimes of individual 
massive stars, clustered star formation produces significantly 
more variability in the number of massive stars present at 
any given time than non-clustered star formation. We de¬ 
fer a discussion of this effect to the following section. Here, 
we instead focus on a second channel by which the CMF 
and fc influence the photometric output, which arises due 
to the method by which slug implements mass constraints. 
As discussed above, a realization of the IMF in a given clus¬ 
ter simulation is the result of the mass constraint imposed 
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Figure 4. Cumulative distribution functions for ionizing luminos¬ 
ity (top panel) and in the GALEX FUV and HST UVIS F225W 
and F336W filters (bottom three panels). In each panel, the x 
axis shows either the ionizing luminosity (for the top panel) or 
the absolute AB magnitude (bottom 3 panels) of a 500 Mq simple 
stellar population drawn from a Chabrier (2005, green), Kroupa 
(2002, blue), or Weidner & Kroupa (2006, red) IMF at ages of I 
Myr (solid), 4 Myr (dashed), and 10 Myr (dot-dashed). 


by the mass of the cluster, drawn from a CMF, and by the 
sampling technique chosen to approximate the target mass. 
Within “galaxy” simulations, a second mass constraint is 
imposed on the simulations as the SFH provides an implicit 
target mass for the galaxy, which in practice constrains each 
realization of the CMF. As a consequence, all the previous 
discussion on how the choice of stopping criterion affect the 
slug outputs in cluster simulations applies also to the way 
with which slug approximates the galaxy mass by drawing 
from the CMF. Through the fc parameter, the user has con¬ 
trol on which level of this hierarchy contributes to the final 
simulation. In the limit of /c = 0, slug removes the interme¬ 
diate cluster container from the simulations: a galaxy mass 
is approximated by stars drawn from the IMF, without the 
constraints imposed by the CMF. Conversely, in the limit of 
/c = 1, the input SFH constrains the shape of each realiza¬ 
tion of the CMF, which in turn shapes the mass spectrum 
of stars within the final outputs. As already noted in Fuma- 
galli et al. (2011), this combined effect resembles in spirit 
the concept behind the IGIMF theory. However, the slug 
implementation is fundamentally different from the IGIMF, 
as our code does not require any a priori modification of the 
functional form of the input IMF and CMF, and each real¬ 
ization of these PDFs is only a result of random sampling 
of invariant PDFs. 

To offer an example that better illustrates these con¬ 
cepts, we perform four “galaxy” simulations with 5000 re- 







Maximum Stellar Mass (Mq) 


Figure 5. Histograms of the actual galaxy mass (top), total stel¬ 
lar mass in clusters (middle), and mass of the most massive star 
formed in clusters (bottom) for four different slug simulations 
with 5000 realizations each. The three simulations labeled as fc 
consist of the default slug CMF and different choices of the frac¬ 
tion of stars formed within clusters. The bottom panel compares 
instead two simulations with /c = 1, but for two different choices 
of CMF (slug’s default and a CMF truncated at 100 Mq). 


alizations. This number ensures that ~ 10 simulations are 
present in each mass bin considered in our analysis, thus pro¬ 
viding a well converged distribution. Each simulation follows 
a single timestep of 2 x 10® yr and assumes a Chabrier (2005) 
IMF, the default ST0P_NEAREST condition, and an input SFR 
of 0.001 M© yr“^. Three of the four simulations further as¬ 
sume a default CMF of the form dN/dM oc between 

20 — 10^ M©, but three different choices of CMF {fc = 1.0, 
0.5, and 0.0). The fourth simulation still assumes /c = 1 
and a dN/dM oc M~^ CMF, but within the mass interval 
20 — 100 Mq (hereafter the truncated CMF). Results from 
these simulations are shown in Figure 5. 

By imposing an input SFR of 0.001 Mq yr“^ for a time 
2 X 10® yr, we are in essence requesting that slug populates 
galaxies with 2000 Mq worth of stars. However, similarly to 
the case of cluster simulations, slug does not recover the in¬ 
put galaxy mass exactly, but it finds the best approximation 
based on the chosen stop criteria, the input CMF, and the 
choice of /c. This can be seen in the top panel of Figure 5, 
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where we show the histograms of actual masses from the four 
simulations under consideration. For /c = 0, slug is approx¬ 
imating the target galaxy mass simply by drawing stars from 
the IMF. As the typical stellar mass is much less than the 
desired galaxy mass, the mass constraint is not particularly 
relevant in these simulations and slug reproduces the tar¬ 
get galaxy mass quite well, as one can see from the narrow 
distribution centered around 2000 Mq. When fc > 0, how¬ 
ever, slug tries to approximate the target mass by means of 
much bigger building blocks, the clusters, thus increasing the 
spread of the actual mass distributions as seen for instance 
in the fc = 0.5 case. For /c > 0, clusters as massive as 10^ 
Mq are potentially accessible during draws and, as a conse¬ 
quence of the ST0P_NEAREST condition, one can notice a wide 
mass distribution together with a non-negligible tail at very 
low (including zero) actual galaxy masses. Including a 10^ 
M© cluster to approximate a 2000 M© galaxy would in fact 
constitute a larger error than leaving the galaxy empty! The 
importance of this effect obviously depends on how massive 
is the galaxy compared to the upper limit of the CMF. In our 
example, the choice of an unphysically small “galaxy” with 
2000 M© is indeed meant to exacerbate the relevance of this 
mass constraint. By inspecting Figure 5, it is also evident 
that the resulting galaxy mass distributions are asymmetric, 
with a longer tail to lower masses. This is due to the fact 
that slug favours small building blocks over massive ones 
when drawing from a power-law CMF. This means that a 
draw of a cluster with mass comparable to the galaxy mass 
will likely occur after a few lower mass clusters have been in¬ 
cluded in the stellar population. Therefore, massive clusters 
are more likely to be excluded that retained in a simulation. 
For this reason, the galaxy target mass represents a limit 
inferior for the mean actual mass in each distribution. 

The middle panel of Figure 5 shows the total amount 
of mass formed in clusters after one timestep. As expected, 
the fraction of mass in clusters versus field stars scales pro¬ 
portionally to /c, retaining the similar degree of dispersion 
noted for the total galaxy mass. Finally, by comparing the 
results of the /c = 1 simulations with the default CMF and 
the truncated CMF in all the three panels of Figure 5, one 
can appreciate the subtle difference that the choice of fc 
and CMF have on the output. Even in the case of /c = 1, 
the truncated CMF still recovers the desired galaxy mass 
with high-precision. Obviously, this is an effect of the ex¬ 
treme choice made for the cluster mass interval, here be¬ 
tween 20 — 100 M©. In this case, for the purpose of con¬ 
strained sampling, the CMF becomes indistinguishable from 
the IMF, and the simulations of truncated CMF and /c = 0 
both recover the target galaxy mass with high accuracy. 
However, the simulations with truncated CMF still impose 
a constraint on the IMF, as shown in the bottom panel. In 
the case of truncated CMF, only clusters up to 100 M© stars 
are formed, thus reducing the probability of drawing stars 
as massive as 120 M© from the IMF. 

This example highlights how the fc parameter and the 
CMF need to be chosen with care based on the problem that 
one wishes to simulate, as they regulate in a non-trivial way 
the scatter and the shape of the photometry distributions re¬ 
covered by slug. In passing, we also note that slug’s ability 
to handle very general forms for the CMF and IMF makes 
our code suitable to explore a wide range of models in which 
the galaxy SFR, CMF and IMF depend on each other (e.g. 



Timestep (Myr) 

Figure 6. Realizations of the SFH in 5000 slug simulations with 
fc = 1. The input SFH is shown in blue, while the dashed black 
line and shaded regions show the median, and the first and third 
quartiles of the distribution. 

as in the IGIMF; Weidner & Kroupa 2006; Weidner et al. 

2010 ). 

3.4 Realizations of a Given Star Formation 
History 

The study of SFHs in stochastic regimes is receiving much 
attention in the recent literature, both in the nearby and 
high-redshift universe (e.g. Kelson 2014; Dominguez Sanchez 
et al. 2014; Boquien et al. 2014). As it can can handle ar¬ 
bitrary SFH in input, slug is suitable for the analysis of 
stochastic effects on galaxy SFR. 

In previous papers, and particularly in da Silva et al. 
(2012) and da Silva et al. (2014), we have highlighted the 
conceptual difference between the input SFH and the out¬ 
puts that are recovered from slug simulations. The reason 
for such a difference should now be clear from the above 
discussion: slug approximates an input SFH by means of 
discrete units, either in the form of clusters (for fc = 1), 
stars (for /c = 0), or a combination of both (for 0 < fc < 1). 
Thus, any smooth input function for the SFH (including a 
constant SFR) is approximated by slug as a series of bursts, 
that can described conceptually as the individual draws from 
the IMF or CMF. The effective SFH that slug creates in 
output is therefore an irregular function, which is the result 
of a superimposition of these multiple bursts. A critical in¬ 
gredient is the typical time delays with which these bursts 
are combined, a quantity that is implicitly set by the SFH 
evaluated in each timestep and by the typical mass of the 
building blocks used to assemble the simulated galaxies. 

A simple example, which also highlights slug’s flexi¬ 
bility in handling arbitrary SFHs in input, is presented in 
Figure 6. For this calculation, we run 5000 slug models with 
default parameters and /c = 1. The input SFH is defined by 
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Figure 7. Whisker plots of the recovered SFH from two “galaxy” simlations of 5000 trials each and input constant SFR of 10“"^ Mq 
yr“^. The cluster fraction was set to fc = 1.0 (left panel) and fc = 0.0 (right panel). At each timestep, the first and third quartiles are 
shown by the box plots, with the median marked by the red lines. The 5— and 95—percentiles are also shown by the whiskers. 


three segments of constant SFR across three time intervals of 
1 Myr each, plus a fourth segment of exponentially decaying 
SFR with timescale 0.5 Myr. Figure 6 shows how the desired 
SFH is recovered by slug on average, but individual models 
show a substantial scatter about the mean, especially at low 
SFRs. An extensive discussion of this result is provided in 
section 3.2 of da Silva et al. (2012). 

Briefly, at the limit of many bursts and small time de¬ 
lays (i.e. for high SFRs and/or when slug populates galaxies 
mostly with stars for /c ~ 0), the output SFHs are reason¬ 
able approximations of the input SFH. Conversely, for small 
sets of bursts and for long time delays (i.e. for low SFRs 
and/or when slug populates galaxies mostly with massive 
clusters for /c ~ 1), the output SFHs are only a coarse repre¬ 
sentation of the desired input SFH. This behaviour is further 
illustrated by Figure 7, in which we show the statistics of 
the SFHs of 5000 galaxy simulations. These simulations are 
performed assuming a constant SFR and two choices of frac¬ 
tion of stars formed in clusters, /c = 1 and 0. One can notice 
that, in both cases, a “flickering” SFH is recovered, but that 
a much greater scatter is evident for the /c ~ 1 case when 
clusters are used to assemble the simulated galaxies. 

From this discussion, it clearly emerges that each slug 
simulation will have an intrinsically bursty SFH, regardless 
to the user-set input, as already pointed out in Fumagalli 
et al. (2011) and da Silva et al. (2012). It is noteworthy 
that this fundamental characteristic associated to the dis¬ 
creteness with which star formation occurs in galaxies has 
also been highlighted by recent analytic and simulation work 
(e.g. Kelson 2014; Dominguez Sanchez et al. 2014; Boquien 
et al. 2014; Rodriguez Espinosa et al. 2014). This effect, and 
the consequences it has on many aspects of galaxy studies 
including the completeness of surveys or the use of SFR in¬ 


dicators, is receiving great attention particularly in the con¬ 
text of studies at high redshift. Slug thus provides a valuable 
tool for further investigation into this problem, particularly 
because our code naturally recovers a level of burstiness im¬ 
posed by random sampling, which does not need to be spec¬ 
ified a-priori as in some of the previous studies. 

3.5 Cluster Disruption 

When performing galaxy simulations, cluster disruption can 
be enabled in slug. In this new release of the code, slug 
handles cluster disruption quite flexibly, as the user can now 
specify the cluster lifetime function (CLF), which is a PDF 
from which the lifetime of each cluster is drawn. This imple¬ 
mentation generalises the original cluster disruption method 
described in da Silva et al. (2012) to handle the wide range 
of lifetimes observed in nearby galaxies (A. Adamo et ah, 
2015, submitted). We note however that the slug default 
CLF still follows a power law of index —1.9 between 1 Myr 
and 1 Gyr as in Fall et al. (2009). 

To demonstrate the cluster disruption mechanism, we 
run three simulations of 1000 trials each. These simulations 
follow the evolution of a burst of 1000 M© between 0—1 
Myr with a timestep of 5 x 10^ yr up to a maximum time 
of 10 Myr. All stars are formed in clusters. The three sim¬ 
ulations differ only for the choice of cluster disruption: one 
calculation does not implement any disruption law, while 
the other two assume a CLF in the form of a power law 
with indices —1.9 and —1 between 1 — 1000 Myr. Results 
from these calculations are shown in Figure 8. 

The total mass in clusters, as expected, rises till a max¬ 
imum of 1000 Mq at 1 Myr, at which point it remains con¬ 
stant for the non-disruption case, while it declines accord- 
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Figure 8. Total mass inside clusters (top) and the galaxy bolo- 
metric luminosity (bottom) as a function of time for three slug 
simulations of 1000 trials each. In all cases, 1000 Mq of stars are 
formed in a single burst within 1 Myr from the beginning of the 
calculation. Simulations without cluster disruption are shown in 
red, while simulations with cluster disruption enabled according 
to a power-law CLF of index —1.9 and —1 are shown in blue and 
green. The black and coloured thick lines show the median, first, 
and third quartiles of the distributions. 


Figure 9. Spectra of four “galaxy” simulations with SFR of 0.001 
Mq yr“^ evaluated at time of 2 x 10® yr. Each simulation differs 
for the adopted extinction law or the inclusion of nebular emis¬ 
sion. In red (solid line), the intrinsic stellar spectrum is shown, 
while models with deterministic and probabilistic extinction laws 
are shown respectively in blue (dashed line) and green (dotted 
line). The solid lines show the first, second, and third quartiles of 
5000 realizations. In brown, the range of luminosities (first and 
third quartiles) with the inclusion of nebular emission is shown. 


ing to the input power law in the other two simulations. 
When examining the galaxy bolometric luminosity, one can 
see that the cluster disruption as no effect on the galaxy 
photometry. In this example, all stars are formed in clus¬ 
ters and thus all the previous discussion on the mass con¬ 
straint also applies here. However, after formation, clusters 
and galaxies are passively evolved in slug by computing the 
photometric properties as a function of time. When a clus¬ 
ter is disrupted, slug stops tagging it as a “cluster” object, 
but it still follows the contribution that these stars make 
to the integrated “galaxy” properties. Clearly, more com¬ 
plex examples in which star formation proceeds both in the 
field and in clusters following an input SFH while cluster 
disruption is enabled would exhibit photometric properties 
that are set by the passive evolution of previously-formed 
stars and by the zero-age main sequence properties of the 
newly formed stellar populations, each with its own mass 
constraint. 


3.6 Dust Extinction and Nebular Emission 

In this section we offer an example of simulations which 
implement dust extinction and nebular emission in post¬ 
processing, two new features offered starting from this re¬ 
lease of the code (see Section 2.1.2). Figure 9 shows the stel¬ 
lar spectra of three slug simulations (in green, blue, and red 
respectively) of a galaxy that is forming stars with a SFR 
of 0.001 Mq yr“^. During these simulations, each of 5000 
trials, the cluster fraction is set to /c = 1 and photometry 


is evaluated at 2 x 10® yr. Three choices of extinction law 
are adopted: the first simulation has no extinction, while the 
second and third calculations implement the Calzetti et al. 
(2000) extinction law. In one case, a deterministic uniform 
screen with Ay = 1 mag is applied to the emergent spec¬ 
trum, while in the other case the value of Ay is drawn for 
each model from a lognormal distribution with mean 1.0 
mag and dispersion of ~ 0.3 dex (the slug default choice). 

As expected, the simulations with a deterministic uni¬ 
form dust screen closely match the results of the non-dusty 
case, with a simple wavelength-dependent shift in logarith¬ 
mic space. For the probabilistic dust model, on the other 
hand, the simulation results are qualitatively similar to the 
non-dusty case, but display a much greater scatter due to 
the varying degree of extinction that is applied in each trial. 
This probabilistic dust implementation allows one to more 
closely mimic the case of non-uniform dust extinction, in 
which different line of sights may be subject to a different 
degree of obscuration. One can also see how spectra with 
dust extinction are computed only for A > 912 A. This is 
not a physical effect, but it is a mere consequence of the 
wavelength range for which the extinction curves have been 
specified. 

Finally, we also show in Figure 9 a fourth simulation 
(brown colour), which is computed including nebular emis¬ 
sion with slug’s default parameters: 0 = 0.73 and log U = 
—3 (see Appendix B). In this case, the cutoff visible at the 
ionisation edge is physical, as slug reprocesses the ionising 
radiation into lower frequency nebular emission according to 


© 2015 RAS, MNRAS 000 , 1-24 



























SLUG III 


15 


10 


7 10 

•< 

7 

to 

P 10 

O) 



37 



10 ^ 



L (erg s'M 


Figure 10. Results of 100 slug and cloudy simulations for a galaxy with continuous SFR of 0.001 Mq yr“^ and fc = 1. Top panel: the 
median SED derived from the 100 slug calculations is shown in red, while the median of the cloudy SEDs computed in integrated mode 
is shown in blue. Bottom panels: histograms of the line luminosity for three selected transitions computed by cloudy in integrated mode 
for each input slug SED. 


the prescriptions described in Section 2.1.2. One can in fact 
notice, besides the evident contribution from recombination 
lines, an increased luminosity at A > 2500 A that is a conse¬ 
quence of free-free, bound-free, and two-photon continuum 
emission. 

3.7 Coupling cloudy to slug Simulations 

In this section, we demonstrate the capability of the 
cloudy_slug package that inputs the results of slug sim¬ 
ulations into the cloudy radiative transfer code. For this, 
we run 100 realizations of a “galaxy” simulation following a 
single timestep of 2 x 10® yr. The input parameters for the 
simulation are identical to those in Section 3.3. The galaxy 
is forming stars at a rate of 0.001 M© yr“^ with = 1. 
We then pipe the slug output into cloudy to simulate H ii 
regions in integrated mode, following the method discussed 
in Section 2.2. In these calculations, we assume the default 
parameters in cloudy_slug, and in particular a density in 
the surroundings of the H ll regions of 10® cm“®. 

Results are shown in Figure 10. In the top panel, the 
median of the 100 slug SEDs is compared to the median of 
the 100 SEDs returned by cloudy, which adds the contribu¬ 
tion of both the transmitted and the diffuse radiation. As 
expected, the processed cloudy SEDs resemble the input 
slug spectra. Photons at short wavelengths are absorbed 


inside the H ll regions and are converted into UV, optical, 
and IR photons, which are then re-emitted within emission 
lines or in the continuum.® The nebula-processed spectrum 
also shows the effects of dust absorption within the H ll re¬ 
gion and its surrounding shell, which explains why the neb¬ 
ular spectrum is suppressed below the intrinsic stellar one 
at short wavelengths. 

The bottom panel shows the full distribution of the 
line fluxes in three selected transitions: Ha, [O III] A5007 
and [N ll] A6584. These distributions highlight how the 
wavelength-dependent scatter in the input slug SEDs is in¬ 
herited by the reprocessed emission lines, which exhibit a 
varying degree of stochasticity. The long tail of the distribu¬ 
tion at low Ha luminosity is not well-characterized by the 
number of simulations we have run, but this could be im¬ 
proved simply by running a larger set of simulations. We are 
in the process of completing a much more detailed study of 
stochasticity in line emission (T. Rendahl et ah, in prepara¬ 
tion). 


® We do not show the region below 912 A in Figure 10 so that 
we can zoom-in on the features in the optical region. 


© 2015 RAS, MNRAS 000 , 1-24 



















16 Krumholz et al. 


3.8 Bayesian Inference of Star Formation Rates 

To demonstrate the capabilities of SFR_slug, we consider 
the simple example of using a measured ionizing photon flux 
to infer the true SFR. We use the library described above 
and in da Silva et al. (2014), and consider ionizing fluxes 
which correspond to SFRq^hO) = 10“^, 10“^, and 10“^ M© 
yr“^ using an estimate that neglects stochasticity and sim¬ 
ply adopts the ionizing luminosity to SFR conversion appro¬ 
priate for an infinitely-sampled IMF and SFH. We then use 
SFR_slug to compute the true posterior probability distri¬ 
bution on the SFR using these measurements; we do so on 
a grid of 128 points, using photometric errors of 0 and 0.5 
dex, and using two different prior probability distributions: 
one that is flat in log SFR (i.e., dp/dlogSFR ~ constant), 
and one following the Schechter function distribution of 
SFRs reported by Bothwell et al. (2011), dp/dlogSFR oc 
SFR^exp(-SFR/SFR.), where a = -0.51 and SFR, = 9.2 
Mq yr“^ 

Figure 11 shows the posterior PDFs we obtain, which 
we normalised to have unit integral. Consistent with the re¬ 
sults reported by da Silva et al. (2014), at SFRs of ~ 10“^ 
M© yr“^, the main effect of stochasticity is to introduce a 
few tenths of a dex uncertainty into the SFR determination, 
while leaving the peak of the probability distribution cen¬ 
tered close to the value predicted by the point mass estimate. 
For SFRs of 10“^ or 10“^ M© yr“^, the true posterior PDF 
is very broad, so that even with a 0.5 dex uncertainty on the 
photometry, the uncertainty on the true SFR is dominated 
by stochastic effects. Moreover, the peak of the PDF differs 
from the value given by the point mass estimate by more 
than a dex, indicating a systematic bias. These conclusions 
are not new, but we note that the improved computational 
method described in Section 2.3 results in a significant code 
speedup compared to the method presented in da Silva et al. 
(2014). The time required for SFR_slug to compute the full 
posterior PDF for each combination of SFRq(hO)? photomet¬ 
ric error, and prior probability distribution is ~ 0.2 seconds 
on a single core of a laptop (excluding the startup time to 
read the library). Thus this method can easily be used to 
generate posterior probability distributions for large num¬ 
bers of measurement in times short enough for interactive 
use. 


3.9 Bayesian Inference of Star Cluster Properties 

To demonstrate the capabilities of cluster_slug, we re¬ 
analyse the catalog of star clusters in the central regions 
of M83 described by Chandar et al. (2010). These authors 
observed M83 with Wide Field Camera 3 aboard the Hubble 
Spaee Teleseope, and obtained measurements for ~ 650 star 
clusters in the filters F336W, F438W, F555W, F814W, and 
F657N (Ho). They used these measurements to assign each 
cluster a mass and age by comparing the observed photom¬ 
etry to simple stellar population models using Bruzual & 
Chariot (2003) models for a twice-solar metallicity popula¬ 
tion, coupled with a Milky Way extinction law; see Chandar 
et al. (2010) for a full description of their method. This cata¬ 
log has also been re-analyzed by Fouesneau et al. (2012) us¬ 
ing their stochastic version of pegagse. Their method, which 
is based on minimization over a large library of simulated 
clusters, is somewhat different than our kernel density-based 
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Figure 11. Posterior probability distributions for the logarith¬ 
mic SFR based on measurements of the ionizing flux, computed 
with SFR_slug. The quantity plotted on the x axis is the offset 
between the true log SFR and the value that would be predicted 
by the “point mass” estimate where one ignores stochastic effects 
and simply uses the naive conversion factor between ionizing lu¬ 
minosity and SFR. Thus a value centered around zero indicates 
that the point mass estimate returns a reasonable prediction for 
the most likely SFR, while a value offset from zero indicates a sys¬ 
tematic error in the point mass estimate. In the top panel, solid 
curves show the posterior PDF for a flat prior probability distri¬ 
bution and no photometric errors (cr = 0), with the three colors 
corresponding to point mass estimates of log SFRgi-jjO) = —5, 
—3, and —1 based on the ionizing flux. The dashed lines show 
the same central values, but with assumed errors of cr = 0.5 dex 
in the measured ionizing flux. In the bottom panel, we show the 
same quantities, but computed using a Schechter function prior 
distribution rather than a flat one (see main text for details). 


one, but should share a number of similarities - see Foues¬ 
neau & Langon (2010) for more details on their method. We 
can therefore compare our results to theirs as well. 

We downloaded Chandar et al.’s “automatic” catalog 
from MAST^ and used duster_slug to compute posterior 
probability distributions for the mass and age of every clus¬ 
ter for which photometric values were available in all five fil¬ 
ters. We used the photometric errors included in the catalog 
in this analysis, coupled to the default choice of bandwidth 
in cluster_slug, and a prior probability distribution that 
is flat in the logarithm of the age and Ay, while varying 
with mass as p(logM) oc 1/M. We used our library com¬ 
puted for the Padova tracks and a Milky Way extinction 
curve, including nebular emission. The total time required 
for cluster_slug to compute the marginal probability dis¬ 
tributions of mass and age on grids of 128 logarithmically- 
spaced points each was ~ 4 seconds per marginal PDF 


^ https://archive.stsci.edu/ 
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Figure 12. Posterior probability distributions for cluster mass 
(left) and age (right) for four sample clusters taken from the cat¬ 
alog of Chandar et al. (2010), as computed by cluster_slug. In 
each panel, the heavy blue line is the cluster_slug result, and 
the thin vertical dashed line is the best fit obtained by Chandar 
et ah. The ID number printed in each pair of panels is the ID 
number assigned in Chandar et al.’s catalog. 


(~ 6000 seconds for 2 PDFs each on the entire catalog of 
656 clusters), using a single CPU. The computation can be 
parallelized trivially simply by running multiple instances of 
cluster_slug on different parts of the input catalog. 

In Figure 12 we show example posterior probabili¬ 
ties distributions for cluster mass and age as returned by 
cluster_slug, and in Figure 13 we show the median and in¬ 
terquartile ranges for cluster mass and age computed from 
cluster_slug compared to those determined by Chandar 
et al. (2010)®. The points in Figure 13 are color-coded by 
the “photometric distance” between the observed photo¬ 
metric values and the 5th closest matching model in the 
cluster_slug library, where the photometric distance is de¬ 
fined as 




-^filter 


\ Vfiit, 


^ ^ (-A^i^obs 


(29) 


where Mi,obs and Mi^uh are the observed magnitude and the 


® The mass and age estimates plotted are from the most up-to- 
date catalog maintained by R. Chandar (2015, priv. comm.). 


magnitude of the cluster_slug simulation in filter z, and the 
sum is over all 5 filters used in the analysis. 

We can draw a few conclusions from these plots. Ex¬ 
amining Figure 12, we see that cluster_slug in some cases 
identifies a single most likely mass or age with a fairly sharp 
peak, but in other cases identifies multiple distinct possible 
fits, so that the posterior PDF is bimodal. In these cases 
the best fit identified by Chandar et al. usually matches 
one of the peaks found by cluster_slug. The ability to re¬ 
cover bimodal posterior PDFs represents a distinct advan¬ 
tage of cluster_slug’s method, since it directly indicates 
cases where there is a degeneracy in possible models. 

From Figure 13, we see that, with the exception of 
a few catastrophic outliers, the median masses returned 
by cluster_slug are comparable to those found by Chan¬ 
dar et al. The ages show general agreement, but less 
so than the masses. For the ages, disagreements come 
in two forms. First, there is a systematic difference that 
cluster_slug tends to assign younger ages for any cluster 
for which Chandar et al. have assigned an age > 10® yr. For 
many of these clusters the 1 — 1 line is still within the 10th 
- 90th percentile range, but the 50th percentile is well be¬ 
low the 1 — 1 line. The second type of disagreement is more 
catastrophic. We find that there are are a fair number of 
clusters for which Chandar et al. assign ages of ~ 5 Myr, 
while cluster_slug produces ages 1 — 2 dex larger. Con¬ 
versely, for the oldest clusters in Chandar et al.’s catalog, 
cluster_slug tends to assign significantly younger ages. 

The differences in ages probably have a number of 
causes. The fact that our 50th percentile ages tend to be 
lower than those derived by Chandar et al. (2010) even 
when the disagreement is not catastrophic is likely a result 
of the broadness of the posterior PDFs, and the fact that 
we are exploring the full posterior PDF rather than select¬ 
ing a single best fit. For example, in the middle two panels 
of Figure 12, the peak of the posterior PDF is indeed close 
to the value assigned by Chandar et al. However, the PDF 
is also very broad, so that the 50th percentile can be dis¬ 
placed very significantly from the peak. We note that this 
also likely explains why our Figure 13 appears to indicate 
greater disagreement between stochastic and non-stochastic 
models than do Fouesneau et al. (2012) ’s Figures 18 and 
20, which ostensibly make the same comparison. Fouesneau 
et al. identify the peak in the posterior PDF, but do not 
explore its full shape. When the PDF is broad or double- 
peaked, this can be misleading. 

The catastrophic outliers are different in nature, and 
can be broken into two categories. The disagreement at the 
very oldest ages assigned by Chandar et al. (2010) is easy 
to understand: the very oldest clusters in M83 are globular 
clusters with substantially sub-Solar metallicities, while our 
library is for Solar metallicity. Thus our library simply does 
not contain good matches to the colors of these clusters, as 
is apparent from their large photometric distances (see more 
below). The population of clusters for which Chandar et al. 
(2010) assign ~ 5 Myr ages, while the stochastic models as¬ 
sign much older ages, is more interesting. Fouesneau et al. 
(2012) found a similar effect in their analysis, and their ex¬ 
planation of the phenomenon holds for our models as well. 
These clusters have relatively red colors and lack strong Ha 
emission, properties that can be matched almost equally well 
either by a model at an age of ~ 5 Myr (old enough for ion- 
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Figure 13. Comparison of cluster masses (left) and ages (right) for clusters in the Chandar et al. (2010) catalog computed using two 
different methods. For each cluster, the value shown on the x axis is the best-fit value reported by Chandar et al. using their non-stochastic 
model grids. The values shown on the y axis are the median masses and ages computed by cluster_slug; for every 20th cluster we also 
show error bars, which indicate the 10th - 90th percentile range computed from the cluster_slug posterior PDF. The dashed black lines 
indicate the 1 — 1 line, i.e., perfect agreement. Points are colored by the photometric distance (see main text for definition) between the 
observed photometry for that cluster and the 5th closest match in the cluster_slug simulation library. 


ization to have disappeared) with relatively low mass and 
high extinction, or by a model at a substantially older age 
with lower extinction and higher total mass. This is a true 
degeneracy, and the stochastic and non-stochastic models 
tend to break the degeneracy in different directions. The 
ages assigned to these clusters also tend to be quite depen¬ 
dent on the choice of priors (Krumholz et al., 2015, in prepa¬ 
ration) . 

We can also see from Figure 13 that, for the most part 
and without any fine-tuning, the default cluster.slug li¬ 
brary does a very good job of matching the observations. 
As indicated by the color-coding, for most of the observed 
clusters, the cluster.slug library contains at least 5 sim¬ 
ulations that match the observations to a photometric dis¬ 
tance of a few tenths of a magnitude. Quantitatively, 90% 
of the observed clusters have a match in the library that is 
within 0.15 mag, and 95% have a match within 0.20 mag; 
for 5th nearest neighbors, these figures rise to 0.18 and 0.22 
mag. There are, however, some exceptions, and these are 
clusters for which the cluster_slug fit differs most dramat¬ 
ically from the Chandar et al. one. The worst agreement 
is found for the clusters for which Chandar et al. (2010) as¬ 
signs the oldest ages, and, as noted above, this is almost cer¬ 
tainly a metallicity effect. However, even eliminating these 
cases there are ~ 10 clusters for which the closest match 
in the cluster_slug library is at a photometric distance of 
0.3 — 0.6 mag. It is possible that these clusters have extinc¬ 
tions Av > 3 and thus outside the range covered by the 
library, that these are clusters where our failure to make the 
appropriate aperture corrections makes a large difference, 
or that the disagreement has some other cause. These few 
exceptions notwithstanding, this experiment makes it clear 
that, for the vast majority of real star cluster observations, 
cluster_slug can return a full posterior PDF that properly 
accounts for stochasticity and other sources of error such 


as age-mass degeneracies, and can do so at an extremely 
modest computational cost. 


4 SUMMARY AND PROSPECTS 

As telescopes gains in power, observations are able to push 
to ever smaller spatial and mass scales. Such small scales are 
often the most interesting, since they allow us to observe fine 
details of the process by which gas in galaxies is transformed 
into stars. However, taking advantage of these observational 
gains will require the development of a new generation of 
analysis tools that dispense with the simplifying assump¬ 
tions that were acceptable for processing lower-resolution 
data. The goal of this paper is to provide such a next- 
generation analysis tool that will allow us to extend stel¬ 
lar population synthesis methods beyond the regime where 
stellar initial mass functions (IMFs) and star formation his¬ 
tories (SFHs) are well-sampled. The slug code we describe 
here makes it possible to perform full SPS calculations with 
nearly-arbitrary IMFs and SFHs in the realm where nei¬ 
ther distribution is well-sampled; the accompanying suite of 
software tools makes it possible to use libraries of slug sim¬ 
ulations to solve the inverse problem of converting observed 
photometry to physical properties of stellar populations out¬ 
side the well-sampled regime. In addition to providing a gen¬ 
eral software framework for this task, bayesphot, we provide 
software to solve the particular problems of inferring the 
posterior probability distribution for galaxy star formation 
rates (SFR_slug) and star cluster masses, ages, and extinc¬ 
tions (cluster_slug) from unresolved photometry. 

In upcoming work, we will use slug and its capability to 
couple to cloudy (Ferland et al. 2013) to evaluate the effects 
of stochasticity on emission line diagnostics such as the BPT 
diagram (Baldwin et al. 1981), and to analyze the properties 
of star clusters in the new Legacy Extragalactic UV Survey 
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(Calzetti et al. 2015). However, we emphasize that slug, its 
companion tools, and the pre-computed libraries of simu¬ 
lations we have created are open source software, and are 
available for community use on these and other problems. 
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APPENDIX A: IMPLEMENTATION DETAILS 
AND CAPABILITIES 

Here we describe some details of slug’s current implementa¬ 
tion. These capabilities may be expanded in future versions, 
but we include this description here both to demonstrate the 


Table Al. Functional forms for PDF segments in slug. 


Name 

Functional form, f{x) 

Parameters^ 

delta 

TT 

1 

o 

xo^ 

exponential 

^ — x/x^ 

x^ 

lognormal 

x~^ exp[—(lnx/a:o)^/(2s^)] 

XQ, s 

normal 

exp[-(a: - XQy/{2s‘^)] 

XQ, s 

powerlaw 

xP 

P 

schechter 

xPe-xtx* 

p, X^ 


^ In addition to the segment-specific parameters listed, all 
segments also allow the upper and lower cutoffs Xa and as 
free parameters. 

^ For delta segments, we require that Xa = X}j = xq. 

code’s flexibility, and to discuss some subtleties that may be 
of use to those interested in implementing similar codes. 


Al Probability Distribution Functions 


Slug uses a large number of probability distribution func¬ 
tions (PDFs). In particular, the IMF, the cluster mass func¬ 
tion (GMF), and the the star formation history (SFH) are 
all required to be PDFs, and the extinction, output time, 
and (for simulations of simple stellar populations) can op¬ 
tionally be described by PDFs as well. In slug, PDFs, can 
be specified as a sum of an arbitrary number of piecewise 
continuous segments. 


dp 

dx 


nifl(x-,Xl,a,Xl,b) + n 2 f 2 {x-,X 2 ,a,X 2 ,b) H-, 


(Al) 


where the normalizations rii for each segment are free pa¬ 
rameters, as are the parameters Xi^a and Xi^h that denote 
the lower and upper limits for each segment (i.e., the func¬ 
tion fi{x) is non-zero only in the range x G [xi^a,Xi^b])- 
The functions fi{x) can take any of the functional forms 
listed in Table Al, and the software is modular so that ad¬ 
ditional functional forms can be added very easily. In the 
most common cases, the segment limits and normalizations 
will be chosen so that the segments are contiguous with one 
another and the overall PDF continuous, i.e., Xi^a = Xi^b 
and nifi{xi^b) = '^?'z+i/z+i(a:z+i,a)- However, this is not a re¬ 
quired condition, and the limits, normalizations, and num¬ 
ber of segments can be varied arbitrarily. In particular, seg¬ 
ments are allowed to overlap or to be discontinuous (as they 
must in the case of S function segments); thus for exam¬ 
ple one could treat the star formation history of a galaxy 
as a sum of two components, one describing the bulge and 
one describing the disc. Slug ships with the following IMFs 
pre-defined for user convenience: Salpeter (1955), Kroupa 
(2002), Ghabrier (2003), and Ghabrier (2005). 

When drawing a finite total mass, in addition to the 
mass distribution, one must also specify a sampling method 
to handle the fact that one will not be able to hit the target 
mass perfectly when drawing discrete objects. Slug allows 
users to choose a wide range of methods for this purpose, 
which we describe briefly here following the naming conven¬ 
tion used in the code. The majority of these methods are 
described in Haas & Anders (2010). 


• ST0P_NEAREST: draw from the IMF until the total mass 
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of the population exceeds Mtarget • Either keep or exclude the 
final star drawn depending on which choice brings the total 
mass closer to the target value. Unless a different scheme is 
deemed necessary, this is the preferred and default choice of 
slug, as this sampling technique ensures that the stochastic 
solution converges towards the deterministic one at the limit 
of sufficiently large Mtarget- 

• ST0P_BEF0RE: same as STOP_NEAREST, but the final star 
drawn is always excluded. 

• STOP_AFTER: same as STOP_NEAREST, but the final star 
drawn is always kept. 

• ST0P_50: same as S'T'0P_1\FEARES'T', but keep or exclude the 
final star with 50% probability regardless of which choice 
gets closer to the target. 

• NUMBER: draw exactly N = Mtarget/(M) objects, where 
(M) is the expectation value for the mass of an object pro¬ 
duced by a single draw, and the value of N is rounded to 
the nearest integer. Note that this method can be used to 
handle the case of characterizing a population as containing 
a particular number of objects as opposed to a particular 
total mass, simply by choosing Mtarget = N{M). 

• POISSON: draw exactly N objects, where the value of N 
is chosen from a Poisson distribution with expectation value 

{N) = Mtarget/(M) 

• SORTED .SAMPLING^: this method was introduced by Wei- 
dner & Kroupa (2006). In it, one first draws exactly N = 
Mtarget/(M) objccts as ffi the NUMBER method. If the result¬ 
ing total mass Mpop is less than Mtarget, the procedure is 
repeated recursively using a target mass Mtarget — Mpop un¬ 
til Mpop > Mtarget- Finally, one sorts the resulting list of 
objects from least to most massive, and then keeps or re¬ 
moves the final, most massive using a ST0P_NEAREST policy. 

Finally, we note two limitations in our treatment of the 
IMF. First, while slug allows a great deal of flexibility in its 
treatment of PDFs, it requires that the various PDFs that 
appear in the problem (IMF, CMF, etc.) be separable, in 
the sense that the one cannot depend on the other. Thus 
for example it is not presently possible to have an IMF that 
varies systematically over the star formation history of a 
simulation. Second, while the IMF can extend to whatever 
mass is desired, the ability of the code to calculate the light 
output depends on the availability of stellar evolution tracks 
extending up to the requisite mass. The set of tracks avail¬ 
able in the current version of slug (Appendix A2) does not 
extend above 120 M©. 


® This method replaces the IGIMF method implemented in the 
earlier version of slug (da Silva et al. 2012), which was based on 
the Weidner et al. (2010) version of the integrated galactic IMF 
(IGIMF) model. In Weidner et al. (2010) ’s formulation, the up¬ 
per cutoff of the IMF in a star cluster depends explicitly (rather 
than simply due to size-of-sample effects) on the total mass of the 
cluster. This model has been fairly comprehensively excluded by 
observations since the original slug code was developed (Fuma- 
galli et al. 2011; Andrews et al. 2013, 2014), and Weidner et al. 
(2014) advocated dropping that formulation of the IGIMF in fa¬ 
vor of the earlier Weidner & Kroupa (2006) one. See Krumholz 
(2014) for a recent review discussing the issue. 


A2 Tracks and Atmospheres 

The evolutionary tracks used by slug consist of a rectan¬ 
gular grid of models for stars’ present day-mass, luminosity, 
effective temperature, and surface abundances at a series of 
times for a range of initial masses; the times at which the 
data are stored are chosen to represent equivalent evolution¬ 
ary stages for stars of different starting masses, and thus the 
times are not identical from one track to the next. Slug uses 
the same options for evolutionary tracks as starburst99 
(Leitherer et al. 1999; Vazquez V Leitherer 2005; Leitherer 
et al. 2010, 2014), and in general its treatment of tracks and 
atmospheres clones that of starburst99 except for minor 
differences in the numerical schemes used for interpolation 
and numerical integration (see below). In particular, slug 
implements the latest Geneva models for non-rotating and 
rotating stars (Ekstrom et al. 2012; Georgy et al. 2013), as 
well as earlier models from the Geneva and Padova groups 
(Schaller et al. 1992; Meynet et al. 1994; Girardi et al. 2000); 
the latter can also include a treatment of thermally pulsing 
AGB stars from (Vassiliadis V Wood 1993). The Geneva 
models are optimized for young stellar populations and likely 
provide the best available implementation for them, but they 
have a minimum mass of 0.8 M©, and do not include TP- 
AGB stars, so they become increasingly inaccurate at ages 
above ~ 10® yr. The Padova tracks should be valid up to 
the ~ 15 Gyr age of the Universe, but are less accurate than 
the Geneva ones at the earliest ages. See Vazquez V Lei¬ 
therer (2005) for more discussion. Models are available at 
a range of metallicities; at present the latest Geneva tracks 
are available at Z = 0.014 and Z = 0.002, the older Geneva 
tracks are available at Z = 0.001, 0.004, 0.008, 0.020, and 
0.040, while the Padova tracks are available at Z = 0.0004, 
0.004, 0.008, 0.020, and 0.050. 

Slug interpolates on the evolutionary tracks using a 
somewhat higher-order version of the isochrone synthesis 
technique (Gharlot V Bruzual 1991) adopted in most other 
SPS codes. The interpolation procedure is as follows: first, 
slug performs Akima (1970) interpolation in both the mass 
and time directions for all variables stored on the tracks; 
interpolations are done in log-log space. Akima interpola¬ 
tion is a cubic spline method with an added limiter that 
prevents ringing in the presence of outliers; it requires five 
points in ID. For tracks where fewer than five points are 
available, slug reverts to linear interpolation. To generate an 
isochrone, the code interpolates every mass and time track 
to the desired time, and then uses the resulting points to 
generate a new Akima interpolation at the desired time. 

Note that the choice of interpolation method does not 
affect most quantities predicted by SPS methods, but it does 
affect those that are particularly sensitive to discontinuities 
in the stellar evolution sequence. For example, the produc¬ 
tion rate of He^-ionizing photons is particularly sensitive 
to the presence of WR stars, and thus to the interpola¬ 
tion method used to determine the boundary in mass and 
time of the WR region. We have conducted tests compar¬ 
ing slug run with stochasticity turned off to starburst99, 
which uses quadratic interpolation, and find that the spectra 
produced are identical to a few percent at all wavelengths, 
with the exception of wavelengths corresponding to photon 
energies above a few Rydberg at ages of ~ 4 Myr. Those dif¬ 
ferences trace back to differences in determining which stars 
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are WR stars, with slug’s Akima spline giving a slightly dif¬ 
ferent result that starburst99’s quadratic one. For a more 
extensive discussion of interpolation uncertainties in SPS 
models, see Cervino & Luridiana (2005). 

Stellar atmospheres are also treated in the same man¬ 
ner as is in starburst99 (Smith et al. 2002). By default, 
stars classified as Wolf-Rayet stars based on their masses, 
effective temperatures, and surface abundances are handled 
using CMFGEN models (Hillier & Miller 1998), those clas¬ 
sified as O and B stars are handled using WM-Basic models 
(Pauldrach et al. 2001), and the balance are treated using 
Kurucz atmospheres as catalogued by Lejeune et al. (1997) 
(referred to as the BaSeL libraries). Different combinations 
of the BaSeL, Hillier & Miller (1998), and Pauldrach et al. 
( 2001 ) atmospheres are also possible. 


APPENDIX B: MODELING NEBULAR 
EMISSION AND EXTINCTION 

Here we describe slug’s model for computing nebular emis¬ 
sion and dust extinction. 


B1 Nebular Continuum and Hydrogen Lines 


As described in the main text, the nebular emission rate can 
be written as La, neb = 7 neb 0 Q(H°)/o^^^(T). Slug computes 
a^^\T) using the analytic fit provided by Draine (2011, his 
equation 14.6), and adopts a fiducial value of cj) following 
McKee & Williams (1997). However, the user is free to alter 
0. Similarly, the temperature T can either be user-specified 
as a fixed value, or can be set automatically from the tabu¬ 
lated cloudy data (see Appendix B2). 

Slug computes the nebular emission coefficient as 


7 .,ne. = 


n<n' 


+a:He7ff^®^ + a:He7M + E tIuL • (B1) 


The terms appearing in this equation are the helium abun¬ 
dance relative to hydrogen XHe, the H^ and He^ free-free 
emission coefficients 7 ^^^ and the H and He bound- 

free emission coefficients 7 ^f^ and 7 bf^\ the H two-photon 
emission coefficient 72 p\ the effective emission rates for H 
recombination lines corresponding to transitions 

between principal quantum numbers n and n , the energy 
differences E^n' between these two states, and the emission 
coefficient for various metal lines (including He lines) 

We compute each of these quantities as follows. For the 
H and He free-free emission coefficients, 7 ^^^ and 7 ^^^^ we 
use the analytic approximation to the free-free Gaunt fac¬ 
tor given by Draine (2011, his equation 10.8). We obtain 
the corresponding bound-free emission coefficients, 7 ^f ^ and 
7bf , by interpolating on the tables provided by Ercolano & 
Storey (2006). We obtain the effective case B recombination 
rate coefficients, by interpolating on the tables 

provided by Storey & Hummer (1995). In practice, the sum 
includes all transitions for which the upper principal quan¬ 
tum number n ^ 25. We compute hydrogen two-photon 
emission via 


(H) _ he . 07 eff,(B),(H) 

72p - ^^(hl )0^2s J 


+ nu/n2s ,crit 


(B2) 


where is the effective recombination rate to the 

2s state of hydrogen in case B, interpolated from the tables 
of Storey & Hummer (1995), n 2 s,crit is the critical density 
for the 2 s — Is transition, and Pj^ is the hydrogen two-photon 
frequency distribution, computed using the analytic approx¬ 
imation of Nussbaumer & Schmutz (1984). The critical den¬ 
sity in turn is computed as 


_ _ A2S-IS 

^2s,crit — ; 7Z ; 7 , 

Q2s — 2p,p T ft + XHej^2s —2p,e 


(B3) 


where A 2 s-is = 8.23 s“^ is the effective Einstein coefficient 
A for the 2s —Is transition (Draine 2011, section 14.2.4), and 
q 2 s- 2 p,p and q 2 s- 2 p,e are the rate coefficients for hydrogen 
2 s — 2p transitions induced by collisions with free protons 
and free electrons, respectively. We take these values from 
Osterbrock (1989, table 4.10). 


B2 Non-Hydrogen Lines 

Metal line emission is somewhat trickier to include. The 
emission coefficients for metal lines can vary over multiple 
orders of magnitude depending on H ii region properties 
such as the metallicity, density, and ionization parameter, as 
well as the shape of the ionizing spectrum. Eully capturing 
this variation in a tabulated form suitable for fast compu¬ 
tation is not feasible, so we make a number of assumptions 
to reduce the dimensionality of the problem. We consider 
only uniform-density H II regions with an inner wind bub¬ 
ble of negligible size - in the terminology of Yeh & Matzner 
( 2012 ), these are classical Stromgren spheres as opposed to 
wind- or radiation-confined shells. To limit the number of 
possible spectra we consider, we also consider only spectral 
shapes corresponding to populations that fully sample the 
IME. Thus while our fast estimate still captures changes in 
the strength of line emission induced by stochastic fluctua¬ 
tions in the overall ionizing luminosity, it does not capture 
the additional fluctuations that should result from the shape 
of the ionizing spectrum. A paper studying the importance 
of these secondary fluctuations is in preparation. 

With these choices, the properties of the H ll region 
are to good approximation fully characterized by three pa¬ 
rameters: stellar population age, metallicity, and ionization 
parameter. To sample this parameter space, we use cloudy 
(Eerland et al. 2013) to compute the properties of H ll re¬ 
gions on a grid specified as follows: 

• We consider input spectra produced by using slug to 
calculate the spectrum of a 10 ^ Mq mono-age stellar popu¬ 
lation with a Ghabrier (2005) IME and each of our available 
sets of tracks (see Appendix A 2 ), at ages from 0—10 Myr 
at intervals of 0.2 Myr. We also consider the spectrum pro¬ 
duced by continuous star formation at M* = 10 “^ Mq yr“^ 
over a 10 Myr interval. Note that the mass and star forma¬ 
tion rate affect only the normalization of the spectrum, not 
its shape. 

• Eor each set of tracks, we set the H ll region metallicity 
relative to Solar equal to that of the tracks. Eormally, we 
adopt cloudy’s “H ii region” abundances case to set the 
abundance pattern, and then scale the abundances of all 
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gas and grain components by the metallicity of the tracks 
relative to Solar. 

• The ionization parameter, which gives the photon to 
baryon ratio in the H ii region, implicitly specifies the den¬ 
sity. To complete this specification, we must choose where 
to measure the ionization parameter, since it will be high 
near the stars, and will fall off toward the H ll region outer 
edge. To this end, we define our ionization parameter to 
be the volume-averaged value, which we compute in the ap¬ 
proximation whereby the ionizing photon luminosity passing 
through the shell at radius r is (Draine 2011, section 15.3) 


Q(H°,r) = Q(H°) 



3 


(B4) 


Here Q(H°) is the ionizing luminosity of stars at r = 0 and 


f 3Q(H°) 
\4:7va^^^n‘^ 


1/3 


(B5) 


is the classical Stromgren radius, which we compute using 
evaluated at a temperature of 10^ K. With these ap¬ 
proximations, the volume-averaged ionization parameter is 


{U) = 


47rr? 




i-ii 


dr 


81 
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1/3 


(B 6 ) 


Thus a choice of (U), together with the total ionizing lumi¬ 
nosity determined from the slug calculation, implicitly sets 
the density nn that we use in the cloudy calculation. Note 
that, as long as nu is much smaller than the critical density 
of any of the important collisionally-excited lines, the exact 
value of Q(H°) and thus nn changes only the absolute lumi¬ 
nosities of the lines. The ratios of line luminosity to (5(H*^), 
and thus the emission coefficients depend only {U). 

Our grid of cloudy calculations uses log {U) = —3, —2.5, 
and — 2 , which spans the plausible observed range. 

After using cloudy to compute the properties of H ll 
regions for each track set, age, and ionization parameter in 
our grid, we record the emission-weighted mean temperature 
(T) = fn‘liTdVlfn‘li dV, and the ratio Liine/Q(H°) for all 
non-hydrogen lines for which the ratio exceeds erg 

photon”^ at any age; for comparison, this ratio is ~ 
erg photon”^ for bright lines such as Ha. This cut typi¬ 
cally leaves ~ 80 significant lines. These data are stored in 
a tabular form. When a slug simulation is run, the user 
specifies an ionization parameter {U). To compute nebular 
line emission, slug loads the tabulated data for the spec¬ 
ified evolutionary tracks and (U), and for each cluster or 
field star of known age interpolates on the grid of ages to 
determine Liine/Q(H°) and thus 7 ^^^; field stars that are 
not being treated stochastically are handled using the val¬ 
ues of Liine/Q(H°) tabulated for continuous star formation. 
At the user’s option, this procedure can also be used to ob¬ 
tain a temperature (T), which in turn can be used in all 
the calculations of H and He continuum emission, and H 
recombination line emission. 

There is a final technical subtlety in slug’s treatment of 
nebular emission. The wavelength grid found in most broad¬ 


band stellar atmosphere libraries, including those used by 
slug, is too coarse to represent a nebular emission line. This 
leads to potential problems with the representation of such 
lines, and the calculation of photometry in narrowband fil¬ 
ters targeting them (e.g., Ha filters). To handle this issue, 
slug computes nebular spectra on a non-uniform grid in 
which extra wavelength resolution is added around the cen¬ 
ter of each line. The extra grid points make it possible to 
resolve the shape of the line (which we compute by adopting 
a Gaussian line shape with a fiducial width of 20 km s“^), at 
least marginally. The grid resolution is high enough so that 
it is possible to compute numerical integrals on the spec¬ 
trum and recover the correct bolometric luminosity of each 
line to high accuracy, so that photometric outputs including 
the nebular contribution can be computed correctly. 


B3 Dust Extinction 

Slug parametrizes dust extinction via the H-band extinc¬ 
tion, Av- As noted in the main text. Ay can either be a 
constant value, or can be specified as a PDF as described in 
Appendix Al. In the latter case, every star cluster has its 
own extinction, so a range of extinctions are present. Once 
a value of Ay is specified, slug computes the wavelength- 
dependent extinction from a user-specified extinction law. 
The extinction curves that ship with the current version of 
the code are as follows: 

• a Milky Way extinction curve, consisting of optical and 
UV extinctions taken from Fitzpatrick (1999), and IR ex¬ 
tinctions taken from Landini et al. (1984), with the two parts 
combined by D. Calzetti (priv. comm., 2014). 

• a Large Magellanic Cloud extinction curve, taken from 
the same source as the Milky Way curve 

• a Small Magellanic Cloud extinction curve, taken from 
Bouchet et al. (1985) 

• a “starburst” extinction curve, taken from Calzetti 
et al. ( 2000 ). 


APPENDIX C: SOFTWARE NOTES 

Full details regarding the code implementation are included 
in the slug documentation, but we include in this Ap¬ 
pendix some details that are of general interest. First, slug, 
cloudy.slug, and various related software packages are fully 
parallelized for multi-core environments. Second, the slug 
package includes a python helper library, slugpy, that is ca¬ 
pable of reading and manipulating slug outputs, and which 
integrates with cloudy_slug, SFR_slug, and cluster_slug. 
In addition to more mundane data processing tasks, the li¬ 
brary supports ancillary computations such as convolving 
spectra with additional filters and converting data between 
photometric systems. The slugpy library is also fully inte¬ 
grated with all the tools described below. FITS file handling 
capabilities in slugpy are provided through astropy (As- 
tropy Collaboration et al. 2013). Third, the code is highly 
modular so that it is easy to add additional options. Extinc¬ 
tion curves and photometric filters are specified using ex¬ 
tremely simple text file formats, so adding additional filters 
or extinction curves is simply a matter of placing additional 
text files in the relevant directories. Similarly, probability 
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distributions describing IMFs, CMFs, SFHs, etc., are also 
specified via simple text files, so that additional choices can 
be added without needing to touch the source code in any 
way. 

The numerical implementation used in bayesphot re¬ 
quires particular attention, since having a fast implementa¬ 
tion is critical for the code’s utility. Since we have written 
the joint and marginal posterior probability distributions of 
the physical variables in terms of kernel density estimates 
(equations 21 and 24), we can perform numerical evaluation 
using standard fast methods. In bayesphot, numerical eval¬ 
uation proceeds in a number of steps. After reading in the 
library of simulations, we first compute the weights from the 
user-specified prior probability distributions and sampling 
densities (equation 18). We then store the sample points 
and weights in a /c—dimensional (KD) tree structure. The 
bandwidth we choose for the kernel density estimation must 
be chosen appropriately for the input library of models, and 
for the underlying distribution they are modeling. There is 
no completely general procedure for making a “good” choice 
for the bandwidth, so bandwidth selection is generally best 
done by hand. 

Once the bandwidth has been chosen, we can evalu¬ 
ate the joint and marginal posterior PDFs to any desired 
level of accuracy by using the KD tree structure to avoid 
examining parts of the simulation library are not relevant 
for any particular set of observations. As a result, once the 
tree has been constructed, the formal order of the algorithm 
for evaluating either the joint or marginal PDF using a li¬ 
brary of Mib simulations is only logWib, and in practice 
evaluations of the marginal PDF over relatively fine grids of 
points can be accomplished in well under a second on a sin¬ 
gle CPU, even for 5-band photometry and libraries of many 
millions of simulations. In addition to evaluating the joint 
or marginal PDF directly on a grid of sample points, if we 
are interested in the joint PDF of a relatively large number 
of physical variables, it may be desirable to use a Markov 
Chain Monte Carlo (MCMC) method to explore the shape 
of the posterior PDF. Bayesphot includes an interface to the 
MCMC code emcee (Foreman-Mackey et al. 2013), allowing 
transparent use of an MCMC technique as well as direct 
evaluation on a grid. However, if we are interested in the 
marginal PDFs only of one or two variables at time (for ex¬ 
ample the marginal PDF of star cluster mass or star cluster 
age, or their joint distribution), it is almost always faster 
to use equation (24) to evaluate this directly than to resort 
to MCMC. The ability to generate marginal posterior PDFs 
directly represents a significant advantage to our method, 
since this is often the quantity of greatest interest. 
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